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OVERVIEW 



The 1980s are a time of profound challenge to the technological strength 
of the United States in the economic as well as the military spheres, and 
our country's performance in research and development in its engineering 
laboratories will be an important and perhaps determining aspect of our 
success in meeting this challenge. Advanced engineering development is 
now based mainly on scientific computation, which in turn relies on math- 
ematical modeling and laboratory experiments. Together they represent 
one of the areas in which the strength of nations is being tested today. 

Mathematics is essential in the development of theoretical and com- 
putational models for solving the highly complex problems of engineering 
and basic science, which encompass a range of scientific difficulties. On one 
side are questions of computer architecture and the science of algorithms. 
On the other side is the modeling of chemical and physical processes 
by means of mathematical equations. The issues are tied together by 
mathematical theory, which seeks a full understanding of the nonlinear 
phenomena contained in the equations and implements this understanding 
through computational methods. This span of scientific activities forms 
the subject that is known as applied mathematics. 

The principal conclusion of this committee is that computational 
modeling, which is a high-leverage element of our nation's scientific and 
technological effort, requires increased emphasis and support. The conclu- 
sion is documented by an examination of typical application areas, which 
reveals the pervasive difficulties that accompany computation of realis- 
tic problems and leads one to consider both what computers can do and 
what they cannot currently do but might eventually be capable of doing. 
As illustration, we examine several deep theoretical problems, including 
turbulence and combustion. At the frontiers of attack on these prob- 
lems we discover the limitations imposed by our current understanding of 
model formulation and computational capability. We examine modeling 
problems and algorithms adapted both to specific features of the desired 
solution and to the computer architecture. We also examine computer 
architecture and component design and manufacture as a mathematical 
modeling problem. 

The Committee recommends increased support for 

1. Research in computational modeling and applied mathematics, 
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2. Computing facilities dedicated to this area, 

3. Education and manpower development in computational and ap- 
plied mathematics. 

These recommended increases include financial support from govern- 
ment and industry as well as institutional support from universities. 



1. INTRODUCTION 



The extensive use of computers in advanced development work began dur- 
ing World War n, and today computing is a vital component of science, 
engineering, and modern technology. Most advanced technological devel- 
opment, from aircraft design to automobiles to petroleum to satellites now 
follows this pattern of reliance on the computer. Moreover, the needs of 
national defense have posed scientific and engineering design problems as 
difficult as any ever encountered. Numerical computation and applied 
mathematics have played an essential role in dealing with such problems. 
In fact, numerical fluid dynamics was born during the 1940s for the pur- 
pose of assisting in the design of nuclear weapons, combat aircraft, and 
conventional ordnance and is now applied widely by industry. 

Most problems of engineering or scientific interest are too complex 
to be modeled and computed exactly. Instead one considers a series 
of approximate models and computations, each of which illuminates a 
different aspect or idealized portion of the overall problem. When used by 
a skilled engineer or scientist, these mathematical models greatly enhance 
the judgment that goes into design decisions and reduce the amount of 
expensive laboratory and field testing required. These advantages account 
for the widespread use of these models. 

More specifically, mathematical models are used in engineering design 
problems in the following modes: 



1. To provide the first test of a new design idea. Beyond common 
sense and simple hand calculations, the computer model is usually the 
cheapest and fastest test of an idea. This test is applied before deciding 
whether to conduct a series of experiments or to build a prototype. 

2. To reduce the time and cost associated with laboratory and 
field tests. Usually engineering problems contain several critical design 
parameters that will have involved a certain amount of trial and error in 
the search for the optimum choice. The computer is used qualitatively 
(Will an increase hi parameter X improve or degrade performance - or 
performance parameter Y"?) and quantitatively (Which values of design 
parameters Xi,.. .,X n will optimize performance?). 

3. To assist in laboratory or field tests that determine model 
parameters and equations. Often the model parameters are measured 



only indirectly. Thus, a computer model of the laboratory apparatus 
may be needed to extract the desired information from the observed data. 
Usually these models are simpler than the complex engineering models, 
and their defining equations can be solved with greater precision and fewer 
approximations. 

4. To replace laboratory or field tests. Sometimes tests are impossible 
or impracticable. For example, measurement of chemical reaction rates 
at extreme conditions of pressure and temperature is very difficult, and 
accumulated experience through trial and error is not adequate for solving 
the problem of landing the Space Shuttle, for reasons of human safety as 
well as cost. 

5. To improve the education and judgment of engineers and scien- 
tists using the models. The mathematical models and computer solutions 
provide a vast increase in the quality and caliber of the data. Thus while 
the laboratory measurement may produce some overall quantity (e.g., total 
flow in and out), the computer model might yield detailed velocities and 
concentrations at each point of the flow field. Because the equations are 
nonlinear, it is difficult to foresee all the relevant phenomena, much less to 
understand then* relative importance. Thus the model becomes for users 
an experimental tool that allows them to understand a problem at a level 
of detail that cannot be achieved by other means. 



In summary, mathematical and computer models are used because 
they are faster, cheaper, and more effective. 

However, models have limitations. There are limitations in the 
validity of the equations used, in the adequacy of the solution algorithm, 
and in the size and speed of the computer. Also, the cost, accessibility, and 
reliability of computer software and, sometimes, the cost of the computa- 
tion itself can be limiting factors. These limitations in some sense define 
the frontiers of science, but more specifically they define the frontiers of 
applied mathematical science. 

Problems of realistic interest typically involve the study of diverse 
physical phenomena oA many scales of length in fully three-dimensional 
settings. Though essential, experimental science in these contexts is ex- 
pensive and difficult. The design of modern strategic weapons systems 
epitomizes these characteristics. For example, in the design of a TRIDENT 
submarine, the architecture of the vessel, the design of the missile, and 
the design of the nuclear warhead all need to be modeled and integrated, 
which requires many thousands of hours of computer time and stretches 
available computing power and modeling techniques to their limits. 

Generally, design and evaluation of new kinds of defense systems, 



such as remotely controlled or robotic vehicles, will require analysis of 
entirely novel systems, using parallel computer architectures among other 
things. The needs of industry for technological advances are similar and 
are dominated by such basic concerns as pollution, depletion of resources, 
energy conservation, and efficient use of manpower. It seems safe to 
say that defense and industrial needs will continue to lead numerical 
computation and applied mathematics into new and challenging regimes. 

Depending on the complexity of the problem and the magnitude of the 
effort expended, models range from excellent to merely suggestive in their 
quality and usefulness. In all cases, improvement of computer modeling is 
one of the most promising avenues to improved technological performance 
by OUT nation. 

As one of the aims of this report, the Committee wants to show and 
emphasize that in the computational approaches to most of today's press- 
ing and challenging scientific and technological problems the mathematical 
aspects cannot and should not be considered in isolation. There is a unity 
among the various steps of the overall modeling process from the formula- 
tion of the physical problem to the construction of appropriate mathemati- 
cal models, the design of suitable numerical methods, their computational 
implementation, and, last but not least, the validation and interpretation 
of the computed results. In particular, the Committee wants to illustrate 
that the steps are more often than not deeply interconnected and that 
the computational process may indeed be part of the model construction. 
At the same tune, there are problem areas, such as turbulence, where cur- 
rent theoretical research may promise a deeper insight into an important 
physical phenomenon. 

In line with these aims, the report uses a "matrix approach" that 
views the same problem from three different standpoints. In Chapter 2, the 
traditional approach is taken of discussing a number of typical problems 
leading to computational modeling from the viewpoint of the scientific 
or engineering discipline in which they arise. Then in Chapter 3 certain 
of these problems are touched on again, this tune from the viewpoint of 
the computational and mathematical difficulties that arise in connection 
with them. For example, these difficulties may involve large numbers of 
degrees of freedom, different scales of time and length, or singularities of 
various types. In Chapter 4 the viewpoint becomes that of the numerical 
algorithms involved in the computations, such as various discretization 
methods, continuation approaches, and splitting techniques. 

Of necessity, many topics have been left out. The list of applications, 
for example, is by no means complete, and, in fact, entire areas such 
as reactor safety and reactor physics are not mentioned at all. Neither 
did the Committee attempt to address all computational and mathemati- 



cal difficulties nor all variations of numerical algorithms. The choice of 
topics was motivated in part by energy-related considerations, the exper- 
tise and interests of the Committee and its advisers, and the report's broad 
purpose, whose achievement would be hindered by any attempt to be en- 
cyclopedic. For the same reason, the report is not intended to be a tech- 
nical summary, and this is also reflected in the fact that no attempt was 
made to reference the relevant literature. The report is mainly addressed 
to scientifically literate readers who know how to consult the literature 
when necesary. 

The Committee obtained advice and technical support from many 
colleagues across the country and from abroad. A list of names of all 
those who helped in this work is given in Appendix A, and the Committee 
is extremely grateful for all the often extensive documentations, special 
write-ups, and other comments that we received. The report is the result of 
a study begun by the Committee in 1981 on computational mathematical 
modeling and mathematics applied to the physical sciences with particular 
reference to the needs of the Department of Energy (DOE). The prepara- 
tion of the report was supported by the Applied Mathematical Sciences 
Research Program of the Office of Basic Energy Sciences of DOE, and the 
Committee also expresses its thanks and appreciation for this support. 

As stressed in the Overview, the Committee found that improvement 
of the mathematical and computer modeling of scientific problems is an 
important priority for our nation. The challenge is broad, and there are 
no simple remedies for current shortcomings. Accordingly, the Committee 
recommends the following: 



1. Increased research support for computational modeling and applied 
mathematics 

The technological challenges of the coming decades will impose new 
tests of our abilities in computational and applied mathematics, and meet- 
ing the tests will require increased research effort. As illustrated in this 
report, the challenges are typically multidisciplinary in nature with applied 
mathematics and modeling often in a central position. 

Hence, to support research in this area, multidisciplinary teams of 
adequate size to make progress on these complex problems should be 
encouraged; and organizational means should be devised to facilitate their 
establishment, continuity, and success. 



2. Increased support for computing facilities dedicated to computa- 
tional modeling and applied mathematics 

Ready access to modern computing systems is essential. Lack of 
equipment is the critical factor most strongly limiting academic research 
in computational mathematics. There is need both for conveniently usable 
local equipment and for access to large-scale computers. Local facilities 
are necessary for entire problems of modest size and for such tasks as 
code development, interactive debugging, test runs, and graphical analysis. 
Large-scale computing is essential because of the size of many typical 
problems as documented in this report. In this connection, the Committee 
strongly endorses the findings and recommendations in the recent Report 
of the Panel on Large-Scale Computing in Science and Engineering, Peter 
D. Lax, Chairman, National Science Foundation, December 26, 1982. 



S. Increased support for education and manpower development in 
computational and applied mathematics 

Today there are unmet manpower needs in computational and ap- 
plied mathematics, as discussed, for example, in Science and Engineering 
Education for the 80s and Beyond, a National Science Foundation report, 
October 1980. These needs are found in industry, government laboratories, 
and academic institutions. The critical challenges in this area call for a 
focus on quality. Specifically, the complex interdisciplinary nature of the 
problems poses special educational challenges for students and young re- 
searchers, and graduate and postdoctoral fellowship support for participa- 
tion in multidisciplinary teams of the type discussed above would also be 
helpful. 

A valuable aspect of such multidisciplinary education is the interac- 
tion it creates between applied mathematicians and other applied scientists 
in universities, government, and industry. 



2. APPLICATIONS 



2.1 HYDRODYNAMIC SYSTEMS 

Hydrodynamic processes touch nearly every aspect of our lives. In a report 
even many times larger than the present one, we could not possibly discuss 
all hydrodynamics applications. Therefore, we have chosen a few, with 
the hope that they will serve to illustrate the importance of hydrodynamic 
models and to point out some of the problems that occur in developing 
the models. 

A large fraction of the current research in computational modeling 
is done with hydrodynamics applications in mind. The great variety of 
responses that we received to our requests for material for this report 
attest to the diversity of these applications. Among them were applications 
having to do with aircraft and wing design, both at subsonic and supersonic 
speeds; global weather prediction and local phenomena such as tornadoes; 
water waves and ship hull design; piping networks, such as in nuclear 
reactor or power plant design; geologic phenomena, such as glacier flow 
or convection hi the Earth's mantle; biological flows, such as the flow of 
blood in the heart; and chemically reacting flow, such as combustion. 

The general system of equations governing hydrodynamics are called 
the Navier-Stokes equations. They are a statement of mass and momen- 
tum conservation, the momentum equation being a formulation of 
Newton's second law, F = ma. The Navier-Stokes equations were first 
developed in France by Navier in the early 1800s. They represented an 
improvement over the Euler equations that were first derived in 1755, in 
that the Navier-Stokes equations included viscous effects that were absent 
in the Euler equations. However, it was not until 1904, when Prandtl de- 
veloped the boundary-layer approximations, that predictions of practical 
viscous flows could be made. Practical solution of the full Navier-Stokes 
equations had to await the development of modern high-speed computers 
beginning in the 1940s. 

Why are equations that have been known for over 100 years so hard to 
solve? The answer lies largely in their inherent nonlinear characteristics. 
Immediately upon looking at the equations one sees that the convective 
transport terms (the acceleration in F = ma) involve velocity times its 
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gradient. This nonlinearity is always present, and it is responsible for 
the existence of complex phenomena such as shock waves and turbulence. 
In principle the Navier-Stokes equations alone provide a description of 
turbulence; however, one would have to resolve such small length scales in 
their solution that this approach is not of practical importance. Therefore, 
a great many approaches to approximating turbulence effects are being 
pursued. Typically, these models introduce further nonlinearities into the 
system. 

Other problems that arise in the solution of hydrodynamic problems 
are related to the disparate time and length scales that must be resolved. 
In particular, convective transport is characterized by the fluid velocity, 
whereas pressure waves travel at the sound speed-typically orders of mag- 
nitude faster than fluid speeds. At the same tune the effects of diffusive 
processes (e.g., shearing stresses) are felt instantaneously throughout the 
flow. In some cases the fluid can react chemically. In those cases, the reac- 
tion rates display a strong nonlinear dependence on the fluid temperature. 
This introduces yet more characteristic scales into the models. 

Depending on characteristic parameters such as the Reynolds number, 
the solutions to the Navier-Stokes equations either can be smooth and 
steady or they can exhibit regular oscillations. Or they can be completely 
chaotic. It is clear that depending on the regime the appropriate solution 
procedures could be quite different. 

Finally a word about boundary conditions. In some problems the flow 
is enclosed, and hence boundary conditions are applied at the boundaries 
of the enclosure (e.g., an automobile engine cylinder), which can often 
have a complex shape. In other cases, such as an aircraft wing, the flow 
is effectively unbounded, and the boundary condition should be applied at 
"infinity" (see Section 3.4). Some modeling problems arise in approximat- 
ing infinity by some finite boundary. In still other problems, such as 
the flow of blood in a heart, the boundary is both complex in shape and 
deforms depending on the forces exerted on it by the fluid. Modeling such 
a problem is clearly a challenge. 



2.1.1 Wings, Wind Tunnels, and Computers 

The economics of the energy shortage implies that planes will fly at speeds 
close to but less than the speed of sound. At such speeds there is a 
"supersonic bubble" over the wing, where the local velocity of the air 
relative to the whig is greater than the speed of sound. In this case the 
presence of shock waves is typical and undesirable. They are undesirable 
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because the drag can be computed as being proportional to the third power 
of the shock strength. The goal of efficient wing design is to produce 
wing shapes with no shocks or only weak shocks in this transonic region. 
A general mathematical theory shows that shockless wing foils exist for 
given transonic cruising speeds. However, the problem of finding such wing 
shapes is both overdetennined and extremely sensitive to small changes in 
the data, i.e., "ill-posed" (see Section 3.6). The solution to such ill-posed 
problems is still valid from an engineering point of view because operation 
at neighboring off-design conditions will produce only weak shocks and 
small drag, see Fig. 2.1. 

Using computer algorithms created by applied mathematicians, it is 
now possible to solve both the inverse problem (design) and the forward 
problem (of determining the flow field for a given wing shape and velocity) 
with sufficient accuracy that the use of costly wind-tunnel experiments can 
be greatly reduced. This accomplishment is a striking success of recent 
applied mathematics. 

The theoretical areas that have contributed to this study include the 
theory of nonlinear elliptic equations, complex function theory, and mixed 
problems, for which a prototype is the Tricomi equation, 



which is elliptic for y > and hyperbolic for y<Q. The elliptic region 
corresponds to the subsonic region, and the hyperbolic region to the su- 
personic "bubble." In a numerical method for the design problem, an 
analytic continuation makes the elliptic region hyperbolic. The resulting 
equation is solved by the method of characteristics, and then the analytic 
continuation back to real values of space is performed numerically. 



2.1.2 Chaos, Turbulence, and Droplets 

Turbulence produces a boundary layer along the trailing edge of an aircraft 
wing. The boundary layer degrades the wing performance and thus is an 
important part of the design problem. The flow behind the trailing edge of 
a wing contains a vortex sheet, and the roll-up of this vortex sheet produces 
turbulence that constitutes a safety hazard for small aircraft flying in the 
wake of a jumbo jet. The axis of the vortex roll-up is perpendicular to the 
wing, and so the roll-up is intrinsically three dimensional. In the simpler 
case of two-dimensional flow the vortex sheet is a line or curve in the 
plane, emanating from the trailing edge of the wing foil. In regions in 
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which the line is stretching, it is geometrically stable. In regions in which 
it is contracting, it is unstable. Instead of contracting, it forms a spiral 
vortex structure and hence is stable. 

The two-phase flow of water and steam in a cooling pipe, or of oil 
and gas in an oil-reservoir production well, is also a problem in which 
the geometrical instabilities of large-scale fluid motion are important. 
Here an internal movable boundary separates regions of different material 
properties. In some cases (a heavy fluid, e.g., water, falling into a light 
fluid, e.g., air), the material interface is unstable against formation of 
fingers. There is continuation of the nonlinear instabilities leading to 
pinchoff of droplets and a chaotic regime (mist) that can be analyzed on 
various length scales, as discussed below in the case of turbulence. 

During the combustion stroke of an automobile engine, the flame 
is quenched when it reaches the cold cylinder walls, and incompletely 
burned fuel present in the combustion chamber at this time contributes to 
pollution and to a loss of fuel efficiency. The rate of advance of the flame 
from the spark plug to the cylinder walls is governed by the laminar flame 
speed and the rate of turbulent mixing. Of these two effects, the second 
is more important. The turbulent mixing is produced by vortices that 
detach from the turbulent boundary layer at the wall during the intake 
and compression strokes (see Fig. 2.2). Thus, an accurate modeling of 
this problem requires an ability to treat a number of fluid singularities: 
flame fronts, vortices, turbulence, boundary layers, and boundary-layer 
separation. 

The examples above show that singularities in fluid flow may be 
geometrically unstable. When this instability occurs in a regime governed 
by the scale-invariant Euler equations, the phenomenon is repeated on 
all length scales and leads to chaotic solutions. Turbulence, vortex roll- 
up, convection fingering, and droplet formation are examples of this 
phenomenon, which we now discuss from a general point of view. 

The Euler equations of fluid dynamics allow intrinsic singularities, 
namely vortices, boundary and shear discontinuity layers, contact or 
material interface discontinuities, and shock waves (Fig. 2.3 shows a com- 
putation of the stretching of a vortex tube in a periodic inviscid flow). 
Depending on the problem, special discontinuities such as flame fronts or 
chemical reaction fronts (for fluid dynamics with chemistry) may occur. 
Within the singularity, the Euler equations fail to be a correct description 
of nature, and corrections (either parabolic effects or perhaps a more com- 
plicated Euler equation with more state variables) may be required. As an 
example, consider a shear layer (i.e., a jump discontinuity in the tangential 
velocity component). Taking the curl of the Navier-Stokes equations, we 
obtain for 
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= V X v = vorticity 



the equation 



+ (v V)u; (u 
at 



where v is the kinematic viscosity. To understand the significance of this 
equation, we specialize to two dimensions. Then w is a vector in the z 
direction, (w V)i>, and 



is the total, or Lagrangian time derivative, so that the Navier-Stokes 
equation says that vorticity moves by passive transport plus diffusion. 
The extra term, (a; V)v above, induces vortex production as a result of 
the stretching of vortex lines in three dimensions and is important for 
considerations of geometrical stability as discussed below. In summary, 
the diffusion term v&v of the Navier-Stokes equation is responsible for 
the vorticity leaving a boundary or internal shear layer and diffusing into 
the rest of the flow. Without viscosity there is no mechanism for vorticity 
to enter (an initially irrotational) flow region. The Prandtl boundary-layer 




FIGURE 2.2 Stretching of a flame by a vortical structure. Such stretching 
is important for the efficient operation of engines; it enhances burning by 
increasing the area of the flame. [From A. J. Chorin, Flame advection and 
propagation algorithms, J. Comput. Phys. 35, 1-11 (1980),] 
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FIGURE 2.3 Successive stages in the stretching of a vortex tube in an 
inviscid periodic flow. Vortex stretching is the mechanism by which energy 
in a turbulent flow is transferred to ever smaller scales where it is eventually 
dissipated, a, step 10, time = 0.65; b, step 20, tune = 0.88; c, step 30, 
time = 1.04; d, step 40, tune = 1.21. [Prom A. J. Chorin, The evolution 
of a turbulent vortex, Commun. Math. Phys. 83, 526-527 (1982).] 
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equations are a special version of the Navier-Stokes equations (scaled in the 
normal direction, so that diffusion occurs only normal to the boundary). 

Often the fluid singularities are geometrically unstable. They may 
bifurcate in a predictable fashion, developing "rolls" (Couette flow) or 
"cells" (Benard flow), or they may become irregular and highly convoluted 
with a tendency toward chaos, known as turbulence. There is no scientific 
reason to question the validity of the Navier-Stokes equation as a micro- 
scopic description of physics even into the turbulent regime. However, its 
usefulness as a description of large-scale fluid motions in the turbulent 
region can be questioned, and some other description of fluid flow, such as 
a random ensemble of interacting vortices, could be more effective. 

There are three energy or length scales in which quite distinct charac- 
teristic phenomena dominate. The smallest length scale is that on which 
energy dissipation dominates. On this length scale, the Navier-Stokes 
equations are the correct equations. The viscosity is large, causing velocity 
fluctuations to be rapidly smoothed and solutions to be (locally) "laminar." 
Since smoothness of solutions is a local property, we may conjecture that 
all solutions of the Navier-Stokes equations should be smooth for all time. 
This conjecture is the major unsolved problem of the energy dissipation 
range. It is known that solutions with smooth data will remain smooth 
for a short time, and solutions with smooth and small data will remain 
smooth for all time. Both statements exclude turbulent regimes. 

For weak solutions, Leray showed that the set of tunes t for which 
v(x, t) fails to be smooth is a set of measure zero. Recently considerable 
progress has been made in restricting the possible singularities of the 
Navier-Stokes equations. 

For longer length scales, viscous effects do not play a major role 
and the fluid flow can be described by Euler equations. However, this 
simplification gives rise to problems. The problems are not merely techni- 
cal but reflect the intrinsically complex phenomenology of fluid dynamics. 
The Euler equations are scale-invariant. Thus if v = v(x, t} is a solution, so 
is v = v(ax, at). The inertia! range is the set of length scales dominated by 
scale-invariant, universal physics. Whatever phenomena can occur (e.g., 
vortices) will be repeated on all length scales in the inertia! range. The 
inertia! range is limited at the smaller end by viscous dissipation. At the 
larger end, it is limited by the special boundary and initial conditions 
imposed on the flow, which result in special (nonuniversal) flow behavior. 

The inertial region is dominated by scale-invariant behavior. There 
is a flow of energy from large-scale motions to smaller ones (the "energy 
cascade"). This cascade seems plausible on physical grounds as a type of 
third law of thermodynamics but does not have a rigorous mathematical 
status. It can be explained as a consequence of the geometrical instability 
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of vortex lines and shear layers. As these go unstable, they generate 
(smaller) new vortices and vorticity. 

The energy cascade leads to a dimensional analysis of characteristic 
exponents and to the Kolmogoroff 5/3 power law 



for the energy distribution as a function of wave number fc. A discussion of 
the experimental data in connection with the Kolmogoroff theory has the 
vortices, which occur on all length scales in the inertia! range, as filling 
space. Actually, it may be better to assume the contrary: vortices of 
a given size fill only a small part of space. Then the smaller vortices, 
which are driven by the larger ones, will occur only within the region of 
these larger vortices, and in fact only within a small part of this region. 
This is the notion of intermittency. It leads to the idea of a singular set 
for solutions of the Euler equations, which is a Cantor set of fractional 
dimension less than 3. 

Intermittency leads to modifications in the Kolmogoroff exponent 
and to a renonnalization group type picture of turbulence. Numerical 
calculations to determine intermittency and energy cascade exponents 
have been performed. The calculations start by tracking vortex lines in 
a fluid flow and proceed through a sequence or renonnalization group 
length scale transformations to focus on the singular Cantor set within 
the solution. 

In most problems, the inertia! region contains lengths too small to be 
used directly in a fluid calculation. Its importance lies hi its role of fixing 
parameters such as an effective or eddy viscosity, which are then used to 
determine the large-scale motion of the fluid. The inertia! region is not 
particularly well understood from either the theoretical or numerical point 
of view. 

The large-scale fluid motions are produced directly by the initial and 
boundary conditions imposed on the flow. These motions are strongly 
problem dependent. Numerical calculations and experiment are important 
tools in their study as is the detailed analysis of simplified and idealized 
flow configurations. An important theoretical question is the evolution of 
initially unstable flow configurations. This question arises in connection 
with the onset of turbulence and in connection with the energy cascade, 
where large-scale vortices excite and drive small-scale ones. 

In some problems (supercritical turbulence), the instability in an ini- 
tially laminar flow is nonturbulent but arises from the bifurcation of a 
fixed point. Further bifurcations lead to higher-dimensional tori, and a 
general theory explains that generically the flow on the (sufficiently high- 
dimensional) torus has a strange attractor as its limit set and that this 
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strange attractor is chaotic in nature. This picture has been analyzed in 
the context of the Lorentz flow, which is the truncation of the Navier- 
Stokes equations to include only a small number of modes. The strange 
attractors found there have a Cantor-like structure. An example of super- 
critical turbulence is Couette flow. 

Subcritical turbulence occurs when the finite (noninfinitesimal) per- 
turbation is less stable than the infinitesimal one. Then turbulence occurs 
below the critical Reynolds numbers at which the linear theory shows 
instability and may go directly to turbulent behavior without a discrete 
sequence of nonturbulent bifurcations starting with laminar flow. 



2.2 CHEMICAL SYSTEMS AND COMBUSTION 

Prom the invention and manufacture of an enormous range of syn- 
thetic materials (e.g., plastics) to the refining and burning of fossil fuels, 
chemistry and chemical processes affect nearly every phase of our lives. 
Naturally it is important to understand and control, as fully as possible, 
many of these complex chemical processes. For example, we seek to find 
new and better materials, to reduce costs, and to generate energy more 
efficiently and with less pollution. Applied mathematics and computa- 
tional modeling continue to play a valuable role in meeting these goals. 

One of the oldest chemical processes harnessed by man is combustion. 
The successful modeling of combustion provides an extraordinarily rich 
source of challenges for the computational mathematician. Frequently, 
combustion models have to incorporate all the difficulties of complicated 
fluid mechanics coupled with complex chemical kinetics. The challenges 
include developing algorithms to ensure accuracy and to reduce computer 
tune and storage. The modeler also seeks appropriate simplifications that 
take advantage of any special attributes of the governing physics in order 
to gain more efficient computation. Since combustion contains a wide 
range of chemical processes, we use it here as an example to illustrate 
points of mathematical interest in general chemical systems. 

Even within the topic of combustion there is an enormous diversity 
of applications. The first topic that probably comes to mind is the model- 
ing of internal combustion engines, and this is an important application. 
Modeling is a part of ongoing research to design new types of engines (e.g., 
direct injected stratified charge), to improve fuel economy, to utilize alter- 
nate fuels (e.g., alcohols), and to reduce pollutant formation. Similarly, 
research for other combustors, such as gas turbines or power-plant boilers 
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benefit from computational models. Still there are many other important 
combustion problems aside from power generation. For example, the field 
of fire research is devoted to problems such as how fires spread in build- 
ings and the behavior of various fire-retardant materials. An important 
current topic in reactor safety is the characterization of hydrogen-air fires 
such as occurred in the Three Mile Island accident. Another public safety 
question deals with the problems of fire and explosion associated with 
a liquid-natural-gas tanker accident. Problems of burning coal and coal 
gasification are also topics of great current interest. 

Perhaps the simplest chemical process from a physical viewpoint is 
chemical equilibrium. At equilibrium all chemical reactions are assumed 
to have gone to completion, and the species concentrations are such that 
the mixture is in a minimum free-energy state. The mathematical com- 
putation of the chemical equilibrium state is posed as a constrained mini- 
mization problem. In combustion the equilibrium composition corresponds 
to the products of combustion long after the combustion is complete. 
Physically, the next step hi complication comes with the inclusion of finite- 
rate chemical reactions. Here the mathematical problem is one of solving 
systems of stiff ordinary differential equations (see Section 3.2 for a discus- 
sion of stiffness). Models of shock tubes or flow reactors, which are used 
frequently to probe fundamental questions in chemical reaction kinetics, 
fall into this class of problems. The physical situation is complicated fur- 
ther by the inclusion of fluid motion and heat and mass transport. In this 
case the mathematical problem is one of solving systems of parabolic or 
elliptic partial differential equations. 

Consider the internal combustion engine as an example. What are 
the things that we might hope to learn from modeling? Ultimately, we 
hope to influence geometrical considerations such as combustion chamber 
shape and component placement (e.g., valves, spark plugs, fuel injectors). 
We also hope to provide a fundamental understanding, on the molecular 
level, about how fuels are oxidized and how pollutants are formed. With 
such understanding we can suggest ways to alter the combustion process 
to advantage. In the past only power and size were important considera- 
tions, and engine optimization could proceed mostly experimentally. Now, 
however, there are too many parameters to optimize simultaneously, and 
computational modeling is increasingly important. 

Because of limitations hi available computational resources, two tacks 
are taken in engine modeling. One is to consider mostly hydrodynamic 
effects. Here the modeling of boundary shapes and component placement 
is of primary importance (e.g., How should the piston face be shaped, 
and where should the spark plug be placed?). Complex domains and two- 
and three-dimensional effects are important. The models must incorporate 
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FIGURE 2.4 Velocity vectors computed in a direct-injection, stratified- 
charge engine at a position near top dead center. The combustion takes 
place in a swirling environment in a cup-like region machined into the 
piston. Several vortices are seen to develop in the cup. [Prom T. D. Butler, 
L. D. Cloutman, J. K. Dukowicz, and J. D. Ramshaw, Multidimensional 
numerical simulation of reactive flow in internal combustion engines, Prog. 
Energy Combust. 3d. 7, 293 (1981).] 

turbulent hydrodynamic effects and sometimes phase-change processes, 
such as fuel spray injections. The chemistry is usually simplified in these 
models because it is not feasible to consider both complex chemical kinetics 
and hydrodynamics on current computers. Figure 2.4 shows the velocity 
vectors that result from a two-dimensional simulation of a direct injected 
stratified charge engine. This is a new engine concept that is being studied 
in a U.S. Department of Energy-sponsored cooperative program including 
General Motors Research Laboratories, Princeton University, and three 
National Laboratories. 
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In addition to the hydrodynamic aspects of engine combustion, there 
are important unanswered questions about the chemistry. Therefore, the 
second tack is to consider simplified hydrodynamic situations, such as 
laminar flames, and treat the chemical kinetics in great detail. These 
models address issues such as ignition phenomena and pollutant formation. 
Figure 2.5 shows some species profiles computed in an atmospheric pres- 
sure acetylene-oxygen premised flame (acetylene is an important reactant 
in soot formation). This model used 30 chemical species and 103 reactions. 
The results were computed using an adaptive mesh placement strategy, 
and they resolve detailed structure within the flame. Note that the flame 
is very "thin"- its thickness is on the order of one millimeter, while com- 
bustion chamber dimensions are on the order of tens of centimeters. 

An interesting aside is to note that these two approaches to combus- 
tion modeling match corresponding approaches to experimental investiga- 
tion. That is to say, it is usually not possible to measure or compute minor 
species concentrations in complex turbulent flows, whereas it is possible 
to do so in laboratory flames. 
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FIGURE 2.5 Species mole fraction profiles showing the internal structure 
of an atmospheric-pressure stoichiometric acetylene-air flame. [From J. A. 
Miller and R. J. Kee, Sandia National Laboratories.] 
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The inherently disparate time scales of chemical reactions, in other 
words their stiffness, contribute to the numerical difficulty of solving com- 
bustion models (see Section 3.2). The inherently disparate tune and length 
scales for fluid transport, heat transfer, and chemical reaction are respon- 
sible for the presence of steep fronts in the solutions. Also, there are many 
degrees of freedom in the system of equations. The number of governing 
partial differential equations is large because a transport equation must be 
included for each species involved in the chemical reaction set. A system 
of 30 to 50 species, involving 100 or more chemical reactions is typical even 
for fuels as simple as methane. Also, for practical combustors, the model 
must ultimately include complex three-dimensional geometries. Because 
complete models of real combustors are too large for present computers, 
an important challenge is to simplify the models (including the physical 
assumptions) to a tractable level. In addition to these problems, there are 
potential difficulties related to scaling. The temperatures are on the order 
of 10 3 K, while some species can have important effects even when their 
mass fractions are as low as 10~ 10 . Moreover, before the computation, 
the peak mass fractions of the various species are usually known only to 
within several orders of magnitude. 

Many of the challenges of combustion modeling have been met by the 
numerical-analysis community; however, many more await resolution. For 
example, for systems of ordinary differential equations we understand how 
to treat the stiffness that results from the complex chemical kinetics by 
using stable implicit methods. However, for systems of partial differential 
equations the application of these methods leaves many open questions 
about how to treat the linear algebra and how to compute the error 
estimates. Operator splitting methods are important in rendering the 
linear algebra tractable for large problems. Stillness also occurs in low- 
Mach-number flows due to very high velocity, but low-amplitude, pressure 
waves. Usually, we do not care about the details of these waves, but 
they can unreasonably limit the size of the time step in explicit methods. 
Subcycling methods, rather than implicit methods, are sometimes used 
to alleviate this problem. Subcycling is a form of operator splitting (see 
Section 4.7) in which the invicid hydrodynamic equations are solved with 
small time steps, while the viscous parts of the equations and the energy 
and mass transport equations are solved using a much larger time step. 

Another particularly important topic in combustion models concerns 
the adequate resolution of localized behavior such as flame fronts. One 
line of research considers adaptive meshing strategies in which a spatial 
mesh network is adjusted dynamically so as to capture the local behavior 
accurately. Other work considers front-tracking methods, where the flame 
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is treated as a local discontinuity and conservation equations connect both 
sides of the front and predict its movement. 

Combustion models must account also for fluid turbulence, an area 
where even the underlying physics is not well understood. Here, computa- 
tional models such as the random vortex method are proving valuable 
in simulating turbulence effects, especially in the investigation of large- 
scale turbulent eddies, the so-called coherent structures. We classify these 
methods as "problem-dependent methods" because the physics and the 
numerical model depend heavily on each other (see Section 4.9). Recently, 
the random vortex method has been combined with a flame propagation 
algorithm so that combustion events can be modeled. Figure 2.6 is a se- 
quence of computer plots that shows the vortex velocity fields and flame 
fronts as computed from a model of turbulent combustion behind a step. 



2.2.1 Asymptotic Analysis 

We mention here that applied analysis (in concert with computation) is 
valuable for many problems hi combustion. Asymptotic methods, for 
example, can actually take advantage of phenomena such as steep fronts. 
They can thus be used to provide insights and to suggest approximations 
that help to reduce the complexity of the numerical models and thus render 
them more tractable. An important aspect of the asymptotic analysis is 
the possibility of considering the dynamical stability of flames. Such work 
may lead to fundamental understanding of such phenomena as the onset 
of turbulence. 

Asymptotic methods exploit the fact that the overall activation energy 
of the chemical reaction is typically large, a well-known consequence of 
which is that flame fronts are very thin. That is, chemical reaction is 
only important in a thin region where the temperature first approaches 
its burned value. On the unburned side of this region, chemical reaction 
is negligible because the temperature is too low, and on the burned side, 
it is negligible because the reaction has essentially gone to completion, 
depleting the unburaed fuel. Mathematically speaking, this thin chemi- 
,cally reactive region may be thought of as an internal boundary layer, 
separating the unburned and burned gases. In the limit, of asymptotically 
large activation energy, the boundary layer is infinitely thin, and we may 
use asymptotic matching principles to connect, or match, the solutions in- 
side and outside the boundary layer. The result is a flame sheet model in 
which the solutions on either side of the sheet are connected by nonlinear 
jump conditions that depend on local conditions at the front. Though the 
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resulting problem is still nonlinear, it represents a significant simplification 
over the original problem with Arrhenius kinetics and is in fact equivalent 
to the original problem in the asymptotic limit of infinitely large activation 
energy. 

The first studies using large activation energy asymptotics began to 
appear about a dozen years ago and were noteworthy for providing for the 
first time analytical flame speed formulas for a steadily propagating planar 
flame. The most significant results to come out of this approach to date are 
the predictions of instability and bifurcation phenomena, which describe 
cellular and pulsating modes of flame propagation. These more complex 
solutions branch, or bifurcate, from the basic solution (in this case, a 
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FIGURE 2.6 A sequential series of computer plots displaying velocity 
fields and flame fronts in turbulent combustion behind a step at inlet 
Reynolds number of 10000. The combustion is stabilized by a recirculation 
region behind the step. The unburned gas is a mixture of propane and air 
at an equivalence ratio of 1/2. [From A. F. Ghoniem, A. J. Chorin, and 
A. K. Oppenheim, Numerical modelling of turbulent flow in a combustion 
tunnel, P/wl Trans, R. Soc. London. Ser. A 304, 303-325 (1982).] 
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steady, planar flame) as one or more parameters exceed critical values for 
instability of the basic solution. Their significance lies partly hi the fact 
that they represent intermediate modes of propagation of the flame in its 
evolution from a laminar to turbulent form of propagation. Experimental 
observations of these transitional modes in their pure form have recently 
been reported, and the theory thus identifies critical parameter thresholds 
that separate one form of flame propagation from another. A number 
of bifurcation analyses using nonlinear perturbation methods have been 
successful in characterizing these nonsteady and/or nonplanar flames for 
parameter values near the instability threshold. However, numerical ap- 
proaches may be used to great advantage in describing these transitional 
solutions when parameters far exceed their critical values for instability of 
the basic solution. 



2.3 PLASMA PHYSICS 

Controlled nuclear fusion provides a possible long-range energy source. 
As currently conceived, a fusion device must contain a deuterium and 
tritium plasma for a sufficiently long tune for net energy production. The 
confinement may be effected either by magnetic fields or by simple iner- 
tia effects. Inertia! confinement of a pellet, coupled with laser or particle 
beam heating, is a distinct possibility, although the bulk of the world 
fusion program centers on magnetic confinement of a plasma. The un- 
derlying problems in magnetic confinement are the determination of the 
equilibrium, stability, transport, and heating properties of plasmas under 
realistic conditions. Numerical modeling and computation represent major 
tools in this study. In recognition of their importance, the Magnetic Fusion 
Energy Computing Center has been established at Lawrence Livermore 
National Laboratory in Livermore, California. This Center is the third 
largest scientific computing center in the United States, and its sole func- 
tion is to provide computing capability to the scientists and engineers in 
the U.S. fusion community. Fusion plasma physics spans a diverse collec- 
tion of fields, with significant efforts occurring in the fields of Hamiltonian 
particle dynamics, statistical mechanics, kinetic theories, and dissipative 
and nondissipative single and multifluid models. The complexity of the 
problems involved makes computation an integral part of the research and 
development program for fusion. 

Controlled-fusion confinement experiments have indicated that, even 
in grossly stable configurations, fluctuations may play an important role in 
detennining energy and particle transport. Thus, an understanding of the 
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nonlinear behavior of various plasma models is fundamental to describing 
such behavior. For both ideal magnetohydrodynamics, a fluid model, and 
the Vlasov-Maxwell system, a kinetic model, linearized equations have 
been studied extensively. These analyses are appropriate for describing 
small-amplitude deviations from a quiescent equilibrium but omit the 
effects of mode-mode coupling and the onset, properties, and evolution 
of turbulence. The development of a self-consistent model of plasma 
equilibrium with fluctuations, stochastic particle and magnetic-field-line 
behavior, and resulting transport will continue to be an important research 
area. 

Plasma fusion applications present many problems in which the equa- 
tions of Hamiltonian dynamics appear. These systems describe single- 
particle phase-space trajectories in the Vlasov-Maxwell theory of colli- 
sionless plasmas. In addition, the trajectory of a magnetic field line in 
a toroidal system is described by a nonlinear Hamiltonian. Thus, the 
question of the existence and construction of adiabatic invariants, ex- 
plicitly time-dependent or not, is of fundamental importance for these sys- 
tems. Significant progress in this area might be made by combining ideas 
from modern topological dynamics, numerical simulation, and perturba- 
tion analysis. For example, if magnetic field lines are ergodic throughout 
a volume rather than lying on closed invariant surfaces, as is given by 
the Kolmorogoff-Arnold-Moser theory, then, owing to electron streaming, 
the thermal conductivity within this volume will be very fast. Thus, such 
questions as when ergodicity arises and what the properties of Hamiltonian 
dynamics are under ergodic circumstances impact strongly on these ap- 
plications. 

Since computational modeling of the full three-dimensional plasma 
equations is out of the question, current work utilizes a range of different 
compromises, simplifying to different degrees either the physics or the 
geometry to obtain a number of computationally tractable problems, each 
of which illuminates in a distinct fashion, a different aspect of the full 
plasma problem. In spite of these difficulties, the computational approach 
is a major route to progress on the problem of controlled fusion. The 
reason is simple. Experiments are expensive and must be supplemented to 
the maximum extent possible by theory. The theory is highly complex and 
nonlinear and is obtained by a combination of numerical experiments and 
physical intuition. Clearly improved computational methods are one of the 
methods through which progress is achieved in this area. Development of 
more powerful computers will also be required, as it is hard to believe that 
smart algorithms will by themselves suffice. 
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2.4 PARTICLE PHYSICS 

The problem is to find the equations that describe the elementary particles 
of subnuclear physics. This problem is among the most difficult to have 
been considered seriously by science hi this century. In recent years there 
have been significant gains in OUT understanding of the mathematical 
structure of the equations of quantum field theory. An analysis of the 
mathematical existence question for the quantum 4 equation 



shows that in space-time dimension d < 3, the theory exists, whereas in 
d > 5, there is no such theory (as would be given by standard methods) 
except for the trivial case X = 0. The method of proof suggests nonex- 
istence for the physical case d = 4 as well and offers confirmation from 
the side of exact mathematical analysis of ideas advanced by theoretical 
physics. This mathematical confirmation is much more compelling than 
any confirmation yet offered from experimental physics. 

The < 4 model is "near" to the Sine-Gordon equation, which is com- 
pletely integrable. Among recent results are the discoveries of a large 
class of stable, localized (in space), periodic (in time) solutions (analogs 
of the Sine-Gordon breather) and of the existence of breather formation 
resonances in soliton/antisoliton scattering. There are many related chal- 
lenging problems in nonlinear mathematics and dynamical systems theory. 

The problem with the 4 equation in d = 4 dimensions is one of several 
reasons for considering in its place quantum gauge fields. A standard 
approach to the quantization problem reduces it to the existence question 
for a singular non-Gaussian functional integral over an mfinite-dimensional 
space. In the case of gauge fields the gauge potentials form an infinite- 
dimensional afGne space A on which the group Q of gauge equivalence 
acts. The functional integral is only defined over the mfinite-dimensional 
manifold 



This manifold is not flat. The integration over M has been shown by exaxit 
mathematical analysis to require introduction of coordinate patches. Only 
locally, within a single coordinate patch can M AfS be regarded as an 
open subset of Euclidean (Hilbert) space. 

One of the methods proposed for understanding integration over M is 
to understand the critical points A 6 M. This question leads to a study of 
classical solutions of the Yang-Mills equation and a reduced form of this 
equation, called the self-dual Yang-Mills equation. Here we are looking at 
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a specific nonlinear elliptic equation in d = 4 dimensions. A remarkable 
analysis has led to a complete classification of its solutions, using methods 
of geometry, topology, and algebra (fiber bundles, the index theorem, and 
algebraic varieties). 

The development of numerical methods for the study of quantum 
fields is in its infancy, mainly because the problem is to a large extent out 
of range of present-day computing machines. Among the methods used 
on this problem, we mention Monte Carlo integration over a space-time 
discretized version of M = AfQ. The discretized problem is called a lattice 
gauge field, and even after discretization, the dimension of M is too high 
to allow evaluation of integrals by direct quadrature. A second approach is 
to generate series coefficients from cluster expansions of statistical physics 
and to reconstruct the desired function space integral from a Fade* analysis 
of the coefficients. Other attempts have been based on finite elements and 
on the renormalization group. It may be that a collaboration of numerical 
analysts with mathematical or theoretical physicists on this problem would 
be beneficial. 



2.5 CONDENSED-MATTER PHYSICS 



2.5.1 Statistical Physics 

In this subject, the equation of state, transport coefficients (pressure, 
viscosity, and thermal conductivity, for example), and other macroscopic 
properties of matter are related to and derived from intermolecular forces. 
The area has diverse and important applications, ranging from metallurgy 
to polymer chemistry to semiconductors and is an active area of research 
from the points of view of theory, numerical methods, and experiment. 
Here we limit ourselves to the two mathematical aspects of this subject: 
mathematical theory and numerical methods. 

The equations of statistical physics involve a large or infinite number 
of degrees of freedom, and so the mathematical theory of use here is 
analysis over infinite-dimensional spaces. Almost the same mathematical 
structure arises as in the study of quantum fields, the relationship between 
the two being that a quantum field is a continuum limit of a statistical 
physics (crystal or lattice) model. 

In a small class of models, including the two-dimensional Ising model, 
exact solutions are known. A larger but still restricted class of models has 
been analyzed mathematically with respect to qualitative behavior. One 
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FIGURE 2.7 Monte Carlo simulation of a kinetic Ising model compared 
with x-ray and neutron deft-action measurements on an alloy of 60 percent 
gold and 40 percent platinum that has been heated and then quenched to 
60 percent of its critical temperature. The abcissa is a reduced momentum 
transfer, and the ordinate is a reduced scattering intensity. (From M. H. 
Kalos, New York University.) 



issue recently addressed hi this work was the stability of surfaces, as is 
relevant to the problem of crystal growth. Another is the effect of aperiodic 
or random crystal structure on bulk material properties. A striking recent 
development was the mathematical demonstration of the dipole binding 
(Kosterlitz-Thouless) phase transition. 

The numerical methods for statistical physics are basically of three 
types. In molecular-dynamics calculations, one takes a large number of 
particles and follows their motion by integration of the ordinary differential 
equation 



dt 2 l ' " '' 

defined by the mtermolecular forces F*. This method is "exact" in its 
treatment of mtermolecular forces, to the extent that they are known 
and that quantum effects can be neglected, but it is approximate in its 
treatment of statistics, since there are computer-dictated limits on the 
number N of particles that can be included. 

28 



Next we mention the method of series expansions, applied to the 
calculation of the equilibrium distribution dw and the partition Z, which 
typically has the form 



and 



z = 



where V is the intermolecular potential energy. Then thermodynamic 
functions such as pressure emerge as derivatives or logarithmic derivatives 
of Z with respect to the parameters (e.g., in V). If the potential energy 



is a sum of pair potentials U then the identity 
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FIGURE 2.8 Energy of liquid helium as a function of density. The broken 
line and points with errors show results of the Monte Carlo quadrature of 
the many-body nonrelativistic Schrodinger equation. The solid line is a fit 
to the experimental data. (From M. H. Kalos, New York University.) 
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substituted in the definition of du above yields an expansion that is rapidly 
convergent, where 

\e- u -l\l 

i.e., where L/~0. This is the region of noninteraction, and so series ex- 
pansions are especially useful to give weak corrections, for example, real 
gas corrections to an ideal gas. The expansions appear to converge up 
to critical points and with considerable work have been used to extract 
information on the equation of state in that region. 

Monte Carlo methods are used in the quadrature of very large- 
dimensional spaces, such as the determination of du above, as well as in 
the direct simulation of stochastic systems. These methods have become 
an experimental tool of mathematical physics. In studying the qualita- 
tive and quantitative behavior of a highly idealized model such as the 
Ising model in equilibrium or very far from equilibrium, or of lattice gauge 
theories, continuum and lattice models of polymers, and atomic models of 
quantum liquids and solids, one can carry out numerical studies in which 
the high-dimensional (i.e., many-body) character of the problem is not 
distorted. 

In application to the calculation of the equilibrium distribution du), 
the essence of the Monte Carlo method is as follows: Starting from an ar- 
bitrary point in a given ensemble, one modules a single particle position Xi 
"at random" but usually so as to lower the potential energy V, occasionally 
so as to raise V. After enough such elementary steps, convergence to the 
distribution duj is obtained. As described here, the method is extremely 
simple, and complications arise from the necessity to obtain convergence 
in a reasonable tune for realistic problems. 

Perhaps the most significant success has been in the microscopic 
theory of classical fluids, where Monte Carlo modeling has provided the 
"experimental" basis for the accurate expansion in \e~ u 1| mentioned 
above. Another conspicuous success is the extraction of critical exponents 
by joining ideas of the renormalization group to Monte Carlo simulation. 
Figure 2.7 shows Monte Carlo simulation of a kinetic Ising model compared 
with scattering measurements on a real alloy. 

Figure 2.8 shows how Monte Carlo methods applied to the many- 
body nonrelativistic Schrodinger equation give a quantitative account of 
the energy of (real) liquid helium as a function of density. 
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2.6 GEOPHYSICAL APPLICATIONS 

Out of a wide range of geophysical applications we limit ourselves here to 
three modeling problems in connection with the discovery and production 
of petroleum. Mathematically the discovery process is an inverse problem 
(see Section 3.6) that of constructing geologic maps using seismic signals, 
which in turn are generated by vibrations or explosions at the surface of 
the Earth. As the signals penetrate the Earth, they cross layers of differing 
density. This causes reflected signals to be returned to the surface, where 
they are then recorded and analyzed. The raw data are very noisy, owing 
to irrelevant near-surface density fluctuations. The noise is removed by 
averaging signals from neighboring receptors or from neighboring source 
locations. Then multiple reflections must be subtracted, and a compen- 
sation (called normal moveout) must be introduced for effects of nonver- 
tical signal propagation. Reflection from nonhorizontal layers generates 
complicating shear waves as well as pressure waves. The subtraction of 
multiple reflection signals can be based on Fourier analysis in a half-space 
and a Wiener-Hopf factorization. 

The correction for nonvertical signals can be based on a reduced 
Helmholtz equation, (-A + fc 2 + v)u - 0, but in order to focus on sig- 
nals moving in only one direction (either up or down), one takes first a 
square root and then a power series expansion of the square root. This 
process is known as the parabolic approximation and leads to the familiar 
Schrodinger equation. Alternatively, the analysis can be based on ray 
tracing and the solution of ordinary differential equations. The numerical 
problem is to implement these steps efficiently, in view of the large amount 
of data to be analyzed. 

The production problem is to describe the flow of oil, water, chemicals, 
and/or heat in a porous sandstone layer. The equations to be solved 
are a coupled system of nonlinear parabolic equations. Generally the 
equations describe mass conservation of individual species, and typically 
some are nearly elliptic whereas others are nearly hyperbolic. In the 
hyperbolic equations, coherent shock and rarefaction waves describing oil 
banks, absorption fronts, and flame fronts may form (see Sections 3.2 and 
3.3). Depending on a dimensionless number known as the mobility ratio, 
heterogeneities, and geometrical effects (convergence versus divergence), 
the fronts may become unstable with respect to the formation of fingers 
(see Section 2.1). 

Critical issues in this problem are the efficient solution of large sparse 
linear systems arising from space and time discretizations and the con- 
trol of numerical instabilities in the solution methods. The physical in- 
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FIGURE 2.9 Successive time steps in the movement of an oil-water in- 
terface toward a producing well. The view is a vertical cross section, 
with the well located on the left boundary of the computational region. 
The instability is caused by three factors: by a heterogeneous conducting 
channel near the bottom of the reservoir, by the converging cylindrical 
flow pattern near the producing well, and by the better flow properties of 
the less viscous displacing fluid (water). [From J. Glimin, B. Lindquist, 

0. McBryan, and S. Yaniv, Statistical fluid dynamics II: the influence of 
geometry on surface instabilities, in Frontiers in Applied Mathematics, Vol. 

1, R. Ewing, ed., SIAM, Philadelphia, Pa. (to be published).] 

stabilities mentioned above are controlled by the use of heat (to make 
the oil flow more easily) or the use of a heavier pushing fluid (polymer- 
thickened water in place of water or water and CO 2 in place of CC>2). 

The third modeling problem we discuss is the process of lifting the 
oil to the surface of the Earth. Because of the sudden drop in pressure, 
considerable gas will typically come out of solution, and the resulting two- 
phase mixture is in the droplet or mist regime mentioned in Section 2.1. 
In contrast to the case of reservoir flow, even the basic equations for two- 
phase flow in a pipe are not well understood. 
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2.7 METEOROLOGY 

Accurate forecasting of the movement of large-scale weather patterns is 
clearly an important problem. Here we mean the tracking of highs and 
lows on the order of 1000 kilometers in diameter. In addition, we are 
also concerned with the prediction of cloud patterns, precipitations, tem- 
perature, and other elements of the weather in as much detail as possible. 
Global predictions of the atmospheric flow are now made routinely at as 
many as 12 levels in the troposphere and stratosphere. 

By studying dynamic models of the atmosphere, estimates have been 
made of the sensitivity of the atmosphere to small perturbations, often 
quantified in terms of the error doubling tune. If two initial states of the 
atmosphere differ by random variations of 1, then it is found that the 
resultant states will differ by 2 in about 3 days and by about 4 in about 
6 days. Hence, meteorologists have said that the weather is unpredictable 
for periods of more than 2 to 3 weeks. But, the forecasters are still far 
from being able to predict the weather accurately for even 1 week. The 
best results are currently obtained by the European Centre for Medium 
Range Weather Forecasting, whose 4-5 day forecasts are currently superior 
to those of the U.S. National Weather Service. 

The accuracy of short-range predictions is limited by four somewhat 
distinct factors: 



1. Initializations: (a) The accuracy and completeness of observational 
data and (b) the compatibility of these data with the mathematical model, 
which is, of course, a simplified representation of the atmosphere. 

2. The limited number of grid points or spectral components used 
in the truncated, i.e., numerical, model as a result of limited computer 
power. 

3. Missing or severely truncated physical processes such as cloud 
dynamics, proper representation of turbulence, and radiation fields, lead- 
ing to inadequacy of the mathematical model describing the atmosphere. 

4. The inherent finite limits of predictability of certain types of non- 
linear dynamical systems such as the atmosphere. (Compare the previous 
paragraph.) 



The study of problems related to the reduction of the influence of 
these factors is at the heart of most research that is directly related to the 
improvement of weather prediction. 

The sparsity and inaccuracy of observational data [l(a)] are being 
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partially overcome by the use of weather-observing satellites. Concerning 
l(b), considerable research is currently directed toward the analysis and 
filtering of the observed data so as to make them compatible with numeri- 
cal models. Items 2 and 3 are closely related since, in a general way, the 
lack of a complete representation of all physical processes is related to the 
incomplete numerical resolution of the mathematical model. Such physi- 
cal processes include cloud dynamics (and the associated cloud physics), 
boundary-layer dynamics (including detailed features of the terrain and 
local heat, moisture, and dust sources), and detailed radioactive processes, 
In item 4 a new concept, "finite limits to predictability" has emerged 
from the work of E. N. Lorenz, namely the idea that the atmosphere and 
perhaps certain other dynamical systems may have finite limits of predic- 
tability regardless of the accuracy and detail of their initial conditions 
and the accuracy with which they are computed. Lorenz concluded from 
a simplified model that hydrodynamical systems with an energy spectrum 
(in the customary sense) having a power-law exponent greater than 3 
(e.g., -5/3) should be unpredictable in detail after a finite time. He sug- 
gests a time of about 15 days for the atmosphere, whose large-scale energy 
spectrum seems to have an exponent of about 2.8, with smaller scales 
(i.e., less than say, 20 km) having the 5/3 power law associated with 
homogeneous turbulence. This question of predictability is clearly related 
to that of the stability of hydrodynamical systems and to modern theories 
of bifurcation and chaos (see Section 3.5). 

A large component of atmospheric research is concerned with topics 
that have little direct application to forecasting. Their motivation may 
lie in basic scientific inquiry or in applications to other disciplines, e.g., 
radio-wave transmission for satellite communications. They may involve 
the study of fundamental hydrodynamic phenomena or other physical 
processes several steps removed from specific application. In all of these 
studies we find extensive use of analytical and numerical modeling of 
physical systems and processes. In addition, a wide variety of statistical 
approaches are used hi all areas of atmospheric research. 

Let us mention a few more areas of study. In addition to the predic- 
tion of large-scale weather patterns, computers are being used in limited 
regions to study the development and movement of smaller-scale distur- 
bances, such as hurricanes, thunderstorms, and tornadoes. Here adaptive 
grid methods are used (see Section 4.2). It is hoped, for example, that a 
better understanding of the details of these destructive phenomena will 
lead to the possibility of altering their course and development. 

On the longer time scale, where seasonal forecasts of average rainfall 
and temperature are made over large regions of the globe, there are at 
present several potentially useful, but unproven ideas. Here, however 
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some new mathematical and physical insight is needed in order to develop 
a satisfactory "average" system of equations. This is an important open 
problem. 

On the still longer time scale of decades and beyond, determining 
the effect of our mechanized civilization on the environment is a basic 
problem. Significant progress has already been made, through the use 
of simple climate models in conjunction with paleontology, astronomy, 
geology, and volcanology, in efforts to understand the factors that influence 
irregular alternations of ice ages and interglacial periods. These statistical- 
dynamical climate models will be essential for predicting the long-term 
effects of the observed increase in carbon dioxide and similar problems 
that will almost certainly arise as our industrial civilization expands. 



2.8 ASTROPHYSICAL APPLICATIONS 

Astrophysical studies pose a wide range of problems for mathemati- 
cal and computational analysis. The underlying physical theories in- 
clude hydrodynamics, magnetohydrodynamics and plasma theory, radia- 
tive transport theory, atomic and nuclear physics, and general relativity 
theory, and hence various comments made in other applications sections 
of this report apply to this area. 

In recent years considerable progress has been made in theoretical and 
numerical studies of a number of astrophysics! topics, as, for instance, the 
study of stellar interiors, the formation of stars and galaxies, the spiral 
structure of galaxies, the physics of supernovae and the evolution of su- 
pemovae remnants, the formation of the solar system, and the behavior of 
binary star systems, to mention just a few. The computational approaches 
may range from relatively simple simulations to the numerical solution of 
complex systems of partial differential equations. 

As an example of the first type, simulations have been applied to 
the study of star formation in a galaxy. It may be assumed that when a 
massive star becomes a supernova the shock wave emanating from it can 
compress the surrounding interstellar gas creating new stars. If at least 
one new star is also a massive star the phenomenon can repeat, leading 
to a chain reaction in the creation of stars. This is equivalent to a direct 
percolation problem, and, as is typical for such problems, phase transitions 
are involved. This is a nonlinear problem with a complicated structure in 
space and time. Accordingly, analytic techniques are difficult, whereas 
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computational simulations are relatively straightforward and lead to the 
development of realistic model galaxies. 

As a second type of example we mention a problem that is amenable 
to considerable mathematical analysis as well as to computational attack. 
This is the question of the spiral structure of galaxies. The so-called 
density wave theory approaches this question as a dynamical problem in 
the form of the gravitational instability of a galactic dish with respect to 
spiral modes. There are three basic approaches hi the calculation of spiral 
structures of galaxies on the basis of density wave theory that may be 
characterized loosely as the stellar model, particle model, and fluid model. 

The stellar model represents the classical approach in the study of 
stellar systems. In brief, a galaxy has a stellar component and a gaseous 
component. The latter usually has a sufficiently small mass to be negligible 
in first approximation. The basic equations governing the behavior of the 
stellar component are the Boltzmann equation, usually in colli sionless form 
when close encounters between stars are omitted, and the Poisson equation 
for the gravitational potential. In a sense, these equations are much simpler 
than the Vlasov equations of plasma physics, but their numerical solution 
nevertheless poses challenging questions. 

In the particle model, the stellar component is considered as a very 
large but finite number of particles. As in our first example, the computa- 
tional approach then assumes the form of a simulation process. In brief, the 
motions of the individual stars are followed, and their gravitational field 
is calculated by a self-consistent evaluation of the field. This technique 
has provided valuable qualitative information. But the number of stars in 
these models is usually of the order of 10 11 - 1 , and hence, by necessity, 
numerical simulations are extremely limited in accuracy and provide only 
few quantitative data for specified galaxy models. 

In the third approach, the stellar system is considered to be a con- 
tinuous medium. In this setting one may then study the characteristics 
of wave patterns over the galactic disk and then dependence on the mass 
distribution. In particular, the spiral structure of galaxies may be ex- 
plained in terms of spiral wave patterns of some kind. Such stationary 
wave patterns over a field of differential rotation are to be expected; in 
fact, hydrodynamic waves over shear flows are known to exist for some 
time especially in the form of self-excited modes. The general form of 
the resulting theory of spiral galactic structures requires elaborate com- 
putational techniques. At the same tune, asymptotic approaches have 
'provided analytic results that have led to a better understanding of the 
dynamical mechanisms despite their limitations in accuracy and in the 
types of galaxies covered. 

As noted earlier, these are only some examples of the many different 
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types of problems in astrophysics that relate closely with applied and com- 
putational mathematics. Many of these problems involve extremely wide 
ranges in length and time scales. For example, calculations of stellar evolu- 
tions have to range over billions of years, while the dynamics, say, of the 
supernovae phase takes place within milliseconds. Similarly, the size of a 
neutron star differs by substantial orders of magnitude from the size of its 
corresponding relevant gravity field. Thus, in particular, the mathematical 
and computational difficulties, discussed in Section 3.2, are especially ap- 
plicable here. Moreover, the wide range of the underlying physical theories 
leads to many substantially different types of mathematical models, which 
in turn require very different computational techniques. 



2.9 STRUCTURAL MECHANICS 

During the past two decades, the use of computers has transformed 
large parts of solid mechanics into practical tools for a multitude of 
technological developments. Sophisticated computational software is 
employed throughout the nation's industries and research laboratories in 
the analysis and design of structures and mechanical equipment. There 
is a strong interaction between applied mathematics and solid mechanics. 
Mathematical analysis has provided insight into model formulation and 
the development of powerful numerical methods; and, vice versa, novel 
engineering approaches have led to new research areas in applied mathe- 
matics. 

In the case of linear problems there exists now a relatively broad 
experience in computing solutions for a range of problems concerning 
the behavior of solid bodies subjected to specified loads. In general, 
the computed results are considered reasonably reliable, and they have 
been corroborated over a period of time by observation and practical 
experimentation. There appears, however, to be a growing need for the 
development of computable error estimates that can provide a realistic 
check on the solution accuracy. Such a posteriori bounds have been shown 
to be feasible; but their applications to large, realistic problems in solid 
mechanics still requires considerable research and software development. 
The general availability of economical a posteriori estimates would make 
possible the consistent use of adaptive mesh-refinement techniques that 
would reduce the cost of data preparation by the users and make it possible 
to generate near-optimal solutions for a given amount of computational 
expense. 
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The general trend in computational solid mechanics today is toward 
extending the computational methodology to nonlinear problems. Sources 
of nonlinearity in structural problems are (1) geometric nonlinearity due 
to nonlinear strain-displacement relations, (2) material nonlinearity due 
to nonlinear constitutive equations, (3) force nonlinearity due to nonlinear 
stress boundary conditions, and (4) kinematic constraint nonlinearity due 
to nonlinear displacement boundary conditions. The source of nonlinearity 
affects the form of the resulting nonlinear equations and, hence, influences 
the effectiveness of the solution techniques. 

The numerical analysis of all of these nonlinear problems is not yet 
at a satisfactory stage. Many computer programs for such problems exist, 
but the mathematical basis for most of the methods used is insufficiently 
understood, and there is little known about the accuracy of the computed 
results. Moreover, in the case of nonlinear problems, few numerical com- 
putations can be supplemented with sufficient experimental experience. 

The situation is still best understood in the case of finite elasticity. 
Even there the mathematical theory of the underlying equations is incom- 
plete, and the approximation theory for these equations is generally based 
on various simplifying assumptions that may or may not be valid for a 
particular problem. 

The state of the art in elasto-dynamic problems is in even worse 
shape. It is known that multiple solutions may exist and shock waves can 
develop. Moreover, not all solutions are necessarily physically relevant. 
The questions of how to model such phenomena numerically and how to 
determine the physically realistic solutions are as yet largely open. When 
it comes to problems in finite plasticity even less is known. Although there 
has been much progress in this area during recent years, no satisfactory 
and complete mathematical model is available as yet. Especially, there 
are profound mathematical and computational difficulties in modeling 
phase changes, viscous effects, cracks and singularities, the growth of 
cracks under dynamic loading, and the identification and implementation 
of physically reasonable constitutive equations describing these materials. 



2.10 NONDESTRUCTIVE TESTING AND TOMOGRAPHIC 
RECONSTRUCTION 

There is a pervasive need in technology to evaluate quantitatively the 
integrity and the remaining reliable lifetime of components and structures, 
e.g., from bridge girders to high-performance ceramic turbine disks. In 
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the past decade considerable progress has been made. This developing 
technology is called nondestructive evaluation (NDE) to distinguish it 
from older nonquantitative nondestructive testing practices. NDE presents 
challenges to the applied mathematical sciences on many levels. A few of 
the more important and topical problem areas are mentioned here. 

In NDE applications a component is often subjected to some sort of 
penetrating radiation with the aim of deducing information about its in- 
ternal state from a measurement of the radiation field external to the part. 
Examples include the use of ultrasonic radiation, x rays, and neutrons. 
Because of its flexibility, relative cost, and safety, ultrasonic methods are 
often used in NDE applications. The deduction of information about flaws 
from the incident and scattered ultrasonic fields relies on the solution or 
approximate solution of the inverse scattering problem for elastic waves. 
In some regimes it is possible to develop and adapt imaging techniques 
to the ultrasonic setting. In either case information about defects must 
be obtained from band-limited noisy data. A fundamental limitation to 
our current ability to utilize ultrasonic techniques broadly is our limited 
understanding of the elastic inverse scattering problem in either the fre- 
quency or the time domain. 

In other important NDE applications, electrical currents are induced 
in a material. These currents produce fields that vary depending on 
whether a defect, e.g., a crack, is present in the material. The utilization 
of these so-called eddy-current techniques depends on the ability to infer 
information about the defects from the measured fields they produce. This 
again is an inverse problem that is imperfectly understood. 

Passive methods in NDS are also widely used. For example, when a 
pressure vessel or aircraft is in service, a crack once formed may grow and 
propagate. This will be accompanied by the release of acoustic energy. 
This energy can in turn be monitored at selected sites. The problem then 
becomes one of identifying and classifying the sources of these sound pat- 
terns using acoustic emission studies, an increasingly important technique 
in NDE. In mathematical terms the problem is the so-called inverse source 
problem, which is beset with the same type of difficulties that the inverse 
scattering problems possesses. Sparse, noisy data often taken at highly 
nonoptimal locations are the raw information from which source charac- 
teristics must be deduced. 

The techniques of NDE have application in other areas, and much 
can be learned in other applications fields that is valuable to NDE. For 
example, there is a close connection between the need for inverse scattering 
results in geophysics and NDE and acoustic imaging results in NDE and 
biomedical applications. An area of great success in medical applications, 
tomography, can provide useful information in selected NDE applications. 
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The greatest success of tomography has been in medical applications 
such as the CAT (computer assisted tomography) scanner. Unlike ordi- 
nary x-ray technique, which masks important features by superposing on 
a single picture information on planes perpendicular to the beam, com- 
puterized x-ray tomography provides pictures of individual thin slices 
through the body. Several hundred parallel x-ray pencil beams are 
projected in the plane of the slice, and the attenuation of the x rays in 
each beam is measured separately and recorded. The procedure is repeated 
for many different beam directions. An elaborate calculation then permits 
approximate reconstruction of the x-ray attenuation density as a function 
of position within the slice. 

The idealized mathematical problem is the reconstruction of a func- 
tion of two variables from its integrals along lines. This problem, as well 
as its three-dimensional version, was solved by Radon in 1917 and later 
rediscovered in various settings such as probability theory (recovering a 
probability distribution from its marginal distributions) and astronomy 
(determining the velocity distribution of stars from the distribution of 
radial velocities in various directions). Of course, much work was needed 
to adapt the Radon inversion formula to the incomplete information avail- 
able in practice. Various algorithms for the numerical inversion of this 
ill-posed problem have been proposed, with the present trend favoring the 
so-called convolution algorithm on account of its speed and accuracy. Each 
area of application has, however, its own requirements and may need a 
modification of existing reconstruction algorithms or even a custom-made 
one. Some algebraic methods, for instance, can easily incorporate a priori 
information about the object to be reconstructed. 

Recent advances in medical tomography include nuclear magnetic 
resonance (NMR) tomography and positron emission tomography (PET). 
In NMR strong magnetic fields are used to affect the nuclear magnetic 
spin rate of hydrogen atoms. By varying the fields and their direction, the 
plane integrals of the density of hydrogen can be measured and the den- 
sity reconstructed by an algorithm based on the above three-dimensional 
version of Radon's theorem. The technique is now regarded as competi- 
tive with x-ray tomography for many purposes; and, of course, it is not 
ionizing The advantage of PET over CAT is that metabolic processes 
can be followed. A compound such as glucose is made using carbon-11 
atoms which emit positrons. Photons resulting from the annihilation of 
an emitted positron with an electron are detected by a bank of detectors 
that can record coincidences. Recently, algorithms based on probabilistic 
arguments have been proposed for the PET reconstruction problem 



40 



2.11 MATHEMATICAL MODELS IN THE BIOLOGICAL SCIENCES 

The biological manifestations of the physical laws of the universe present 
us with a rich variety of new phenomena that require the devlopment 
of new mathematical tools and computational methods. We shall discuss 
just a few examples of mathematical research in cardiovascular physiology 
and neurophysiology, with the knowledge that there are many other areas 
of biological sciences in which mathematics and computing are fruitfully 
applied. 

Blood flow in the heart obeys the incompressible Navier-Stokes equa- 
tions, which, in turn, are simply a statement of Newton's laws in 
differential form (see Section 2.1). The distinctive biological character 
of the problem comes, however, from the moving boundaries that are in 
contact with the blood. These include the muscular heart walls and the 
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FIGURE 2.10 Computer-generated plots showing the predicted opening 
movement of a curved pivoting-disk valve mounted in the mitral (inflow) 
position of the left ventricle. The curvature makes the valve open more 
widely than a straight valve pivoted at the same point. It also helps 
to prevent stagnation in the smaller opening of the valve. [Prom D. M, 
McQueen and C. S. Peskin, Computer-assisted design of pivoting-disc 
prosthetic mitral valves, J. Thorac. Cardiovaac, Surg. (in press),] 
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1.4 

FIGURE 2.11 Transmission and "reflection" of a nerve impulse at a junc- 
tion where the diameter of the neuron suddenly increases. The plots show 
computed voltage as a function of tune at equally spaced positions. The 
junction is at x = XQ, and the ratio of diameters there is 2.5:1. Note the 
increase in propagation speed for x > XQ. A reflected wave is set up when 
the larger fiber re-excites the smaller fiber after the refractory period of the 
smaller fiber has elapsed. [From S. S. Goldstein and W. Rail, Changes of 
action potential shape and velocity for changing core conductor geometry, 
Biophysical J. 14, 731 (1974).] 



elastic heart valve leaflets. The motions of these boundaries are not known 
in advance; they must be computed along with the motion of the fluid. 

These considerations have led to the development of a computational 
model of the left heart that can be used in the computer-aided design 
of artificial heart valves. In this model, the fluid equations are solved 
by finite-difference methods on a regular, square mesh (see Section 3.4). 
The boundaries are represented in Langrangian form as a collection of 
moving points. Coupling coefficients between boundary markers and fluid 
mesh points are computed with the aid of an approximation of Dirac's 
^-function. This computer model (Fig. 2.10) has been used in the design 
of prosthetic heart valves to remedy problems of stagnation and blood 
clotting in the smaller opening of the valve. The model has also been 
helpful in the study of disease processes, providing, for example, a possible 
explanation for mitral valve prolapse. 
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Just as the flow of blood in the heart is ultimately governed by 
Newton's laws, the conduction of electrical signals along nerves is ul- 
timately governed by Maxwell's equations (Figs. 2.11 and 2.12). Here 
also, nature provides a peculiar boundary condition that leads to entirely 
new phenomena. In this case the boundary condition comes from certain 
voltage-dependent channels located in the nerve membrane. These intro- 
duce a nonlinearity, and the equations of nerve conduction take the form 
of nonlinear diffusion equations the Hodgkin-Huxley equations. Without 




FIGURE 2.12 Computed electrical signals at the input end of a neuron. A 
brief pulse of current is applied at the periphery of a tree, and the result- 
ing voltages are computed (logarithmic scale) at the input (BI), at succes- 
sive branch points (P, GP, GGP), and finally at the cell body (SOMA). 
[Modified from J. Rinzel and W. Rail, Transient response in a dendritic 
neuron model for current injected at one branch, Biophysical J. 14, 759 
(1974).] 
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the nonlinearity, signals introduced at one end of the nerve would decay 
rapidly; because of the nonlinearity such signals evolve into a specific 
waveform that propagates at a constant speed without distortion. This 
waveform is called the nerve impulse, and it is the basic unit of long- 
distance communication in the nervous system. 

Currently, there is a wide range of mathematical, computational, and 
physiological research activity related to the Hodgkin-Huxley equations. 
Mathematically, there is extensive research on the basic theory of non- 
linear diffusion equations. A particularly fruitful approach here has been 
the use of piecewise linear models that expose the basic structure of the 
equations. Singular perturbation methods have also been useful because 
the equations exhibit a disparity of time scales. 

An important physiological enterprise is the modification and applica- 
tion of the Hodgkin-Huxley equations to other excitable tissues. In the 
heart, for example, equations of the Hodgkin-Huxley type describe the 
electrical processes that generate the cardiac rhythm and coordinate the 
heartbeat. Mathematicians are just beginning to use these equations as 
a basis for a theory of the abnormal rhythms of the heart, of which the 
most serious is ventricular fibrilation. This theory has connections with 
recent work on chaotic dynamical systems and transition to turbulence: it 
appears that fibrilation is directly analogous to turbulence. This work has 
enormous practical significance, since the principal cause of sudden death 
following heart attacks is ventricular fibrilation. 

Progress has also been made in the modeling of the input to the 
neuron (whose output signal is the nerve impulse). The neuron integrates 
information received through a tree of dendrites in which the signal- 
ing mechanism is often described by the linear diffiision equation with 
leakage. Mathematical modeling of the dendritic tree has had a substan- 
tial impact on experimental neurophysiology. One reason for this is that 
dendrites are too small to be penetrated with microelectrodes. Thus the 
neurophysiologist can only record voltage or inject current at the cell body 
and is forced to rely on the theory to indicate the significance of these 
measurements with respect to activity in the dendritic tree. 

Some major successes of the theory are as follows: 

1. Elucidation of the dramatic differences between effects of a synapse 
close to the cell body and effects of a similar synapse far out in the dendritic 
tree. 

2. A possible explanation of the role of dendritic spines in learning 
and memory. 

3. Prediction of the existence of dendro-dendritic synapses based 
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on a mathematical model of field potentials in the olfactory bulb. Such 
synapses, previously unheard of, were subsequently found in electron 
micrographs. This work led to a fundamental new concept of local in- 
formation processing in an "analog" rather than a "digital" mode, i.e., 
without nerve impulses. Such processing is important in neural networks 
such as the retina. Indeed, an extensive mathematical theory of the retina 
has been developed, and this is another exciting area of current research. 



2.12 ELECTRONIC COMPONENTS 

The design and fabrication of modern integrated circuits is a complex 
process. The number of devices that one can put on a chip depends on 
the size of the chip and how small one can make its features. Over the 
years the largest increase in the number of devices on a chip has resulted 
from the continuing reduction in feature size and with this a reduction in 
device size. 

Consequently, process and design engineers have had continually to 
redesign the process steps and then recalculate the resulting device charac- 
teristics to ensure good electrical behavior. This has had to be computa- 
tional because the equations are mathematically intractable and a trial and 
error approach is prohibitively expensive and tune-consuming. Moreover, 
experimental techniques tell us only what happened, not why. Effective 
device design depends on determmmg both the what and the why by vary- 
ing the problem parameters in the computational model. 

The mathematical models on which the theory of semiconductor 
devices rests are differential equations that describe the flow of current 
(holes and electrons) under the influence of electric fields. When feature 
sizes were large, the devices could be treated as though they consisted 
of plane surfaces and edge effects could be neglected. This allowed the 
development and successful use of one-dimensional analytical models ob- 
tained by solving a system of three coupled, nonlinear ordinary differential 
equations. As device sizes shrunk, these models became more complicated 
and less accurate as edge effects became more important. In very large- 
scale integration where device dimensions have reached a few micrometers, 
these models are no longer adequate, and the coupled, nonlinear, partial 
differential equations must now be solved in two and sometimes three 
dimensions. These differential equations consist of a nonlinear Poisson 
equation that describes the potential of the electric field and two non- 
linear transport equations that describe the motion of the holes and the 
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electrons. Moreover, the smaller sizes have made some physical effects 
important whereas they were previously neglected. This results in further 
complications to the basic equations. 

Partial differential equations play an increasingly important role 
in simulating the fabrication processes. The transistors in a chip are 
formed by implanting certain dopant ions into selected areas of the chip. 
Subsequent high-temperature processes, such as growing an oxide layer, 
will cause these atomic impurities to diffuse. Their final distribution is an 
important factor in determining device characteristics. Insights into these 
fabrication processes are especially important for increasing the yield of 
reliable devices, which is a critical factor in their economic viability. 

The design of the overall circuit to be placed on a chip leads to large 
systems of nonlinear differential equations that need to be solved numeri- 
cally. Then the efficient layout of the circuit on the chip introduces com- 
binational and graph theoretical problems, which again pose formidable 
computational problems. 
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3. COMPUTATIONAL AND MATHEMATICAL DIFFICULTIES 



Computational mathematics and modeling is a relatively young science, 
one that is expanding rapidly. The success and importance of the field 
stems from the fact that its application provides the possibility of tack- 
ling significantly more complex and difficult problems than would other- 
wise be possible. Perhaps the greatest opportunity provided by a com- 
putational approach is that it opens the wide realm of strongly non- 
linear phenomena to systematic, relatively accurate, and efficient model- 
ing, improving the chance that important phenomena can be isolated 
and analyzed. Nonlinearities pervade nearly all aspects of applied mathe- 
matics, and to a large extent these nonlinearities are responsible for the 
difficulties that are encountered in computational modeling. Our purpose 
in this chapter is to explore the source of some of these difficulties. In 
Chapter 4, we will discuss the direction of some of the computational 
research needed to resolve the problems. 



3.1 DEGREES OF FREEDOM 

There are several reasons for the degrees of freedom in a model, and 
hence the size of computational problems, to increase. One is that we 
attempt to increase the accuracy and complexity of our representations of 
the physical conservation laws. Increased degrees of freedom come either 
from increasing the number of dependent or independent variables. For ex- 
ample, in chemical models the number of dependent variables is increased 
by increasing the number of chemical species considered. An obvious 
need for increased independent variables comes from the need to represent 
phenomena in two and three spatial dimensions. However, even higher- 
dimensional problems arise when the independent variables are not the spa- 
tial coordinates but are various state descriptors; such higher-dimensional 
problems are common in physics and chemistry, Unfortunately, the direct 
application of numerical methods that work well in one or two dimensions 
often are not usable in three dimensions. Therefore, increasing the degrees 
of freedom may require significantly different algorithms. 
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Another way to increase the degrees of freedom in problems is to 
allow the possibility of multivalued solutions. Certain physical systems 
allow the solutions to bifurcate, a phenomenon discussed in Section 3.5. 
Other systems have hysteresis or "history-dependent" properties. In these 
problems the solutions depend not only on the boundary conditions but 
also on the path that the transient solution follows. Structural dynamics 
models for plastic materials in which the stress depends not only on the 
strain but on the straining rate as well exhibit this behavior. 

Consider an example to quantify the magnitude of the computational 
requirements for multidimensional problems. First suppose we want to fol- 
low the evolution of a chemical model having N chemical species. Suppose 
also that the problem is stiff (see Section 3.2) so for each time step in 
the evolving system we solve N nonlinear algebraic equations in N un- 
knowns. .Even for models with a large number of species this is a relatively 
straightforward task. However, if we now want to introduce transport 
phenomena, such as fluid mixing, the problem has to include the spatial 
dependence of each chemical species. 

For a one-dimensional case assume that we use I spatial mesh points 
and that we estimate what each species will do on the basis of its current 
local value and that of its immediate neighbors (a three-point spatial stencil 
involving all species). We now have NI unknowns and a nonlinear system 
that has 3NI nonzero entries in each equation. It is typical to have 40 to 50 
species and to require 100 mesh points to resolve the species concentrations 
accurately. In this case over three quarters of a million words of memory 
are needed just to store the approximating local linear system along with 
the solution. This alone is larger than the fast-access memory of most 
modern computers. Suppose further that the same model is to be posed 
in two dimensions on a I x I mesh and in three dimensions on an / X / X / 
mesh. The solutions themselves require NI 2 and NI 3 words of storage, 
respectively. Worse, the totality of coefficients in the typical linear system 
will be 5N 2 / 2 and 7N 2 I 3 . It is of practical interest to want solutions to 
three-dimensional combined kinetics and transport problems. However, 
even for a modest system of 20 species and 50 mesh points per spatial 
dimension, 2.5 million words are needed to store the solution alone. An 
additional 350 million are required to store the coefficients of the linear 
system. (For comparison, we note that the largest computers currently 
available have 4 million words of fast-access memory.) Moreover, in the 
high-dimensional cases, the linear system is not conveniently structured 
for efficient solution. 

Clearly a major problem in solving high-dimensional systems of par- 
tial differential equations is that after discretization the resulting system 
of approximating linear equations can be much too large to be solved 

48 



effectively by direct means. Except in simple cases the resulting equations 
have to be solved iteratively, and, especially for strongly nonlinear equa- 
tions, the convergence properties of the iterative process may be a major 
concern. 

Operator splitting methods, such as alternating direction implicit 
(ADI) (see Section 4.7), also provide a current effective approach for 
high-dimensional problems. The approach here is to alternate the solu- 
tion of a set of lower-dimensional problems. The alternated problem ap- 
proximates the original problem with sufficient accuracy, but the set of 
lower-dimensional problems is much more easily solved, even in aggregate. 
Again convergence and accuracy properties must be established in all but 
the simplest of cases. 

For some problems it may be more efficient to depart from the con- 
ventional ideas of discretization on a mesh network and consider instead 
mesh-free methods or (see Sections 4.8 and 4.9) Monte Carlo methods. 
Unknowns can be changed: instead of the amount of material at a given 
location, one can ask for the amount of a given wavelength in the solution 
as a whole. In fluid mechanics, vortex methods could be a more efficient 
approach. Such methods reduce the size of the linear algebra problems, in 
comparison with the mesh-oriented methods, if the new unknowns carry 
all the important information desired. 

Decisions about which methods are most effective may depend 
strongly on advances and changes in computer architecture. Since new 
architectures, such as vector and parallel processors, are now evolving, 
the numerical analyst has to re-evaluate his approaches periodically. Also, 
considerations such as the relative cost of memory versus central processor 
time can bear heavily on decisions regarding algorithms. 



3.2 DIFFERENT TIME AND LENGTH SCALES 

In principle a model can contain important length scales that range from 
the size of an atom to the size of the universe. In practice, however, we 
limit the range of scales by approximation (see, for example, Section 3.7). 
Nevertheless, it is often the case that mathematical models of physical 
processes are characterized by the simultaneous presence in their solution 
of significantly different time and length scales. The solutions to these 
models will have regions of strongly localized behavior, such as shocks, 
steep fronts, or other near discontinuities. Therefore important topics of 
research in numerical analysis are the consideration of such circumstances 
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and the development of efficient theories and methodologies for their com- 
putational solution. Indeed, we often find situations for which solutions 
are not possible, or at least not practical, without the application of spe- 
cialized methods to deal with the multiple characteristic scales. 

Fluid mechanics and chemical kinetics are two areas that provide a 
rich source of examples for multiple and disparate scales. Fluid-mechanical 
processes are commonly characterized by groups of parameters that are 
indicative of the various scales in a problem. Some examples are the Mach 
number (relates velocities to sound speed), the Reynolds number (relates 
inertia! forces to viscous forces), the Prandtl number (relates viscous effects 
to thermal effects), and the Damkb'hler number (relates chemical reaction 
rates to diffusion rates). When any of these numbers is very large (or 
small) it is likely that the solutions to the models will have regions of 
localized behavior. For the case of large Mach number the possibility of 
shocks exists. Similarly for large Reynolds number we expect to encounter 
boundary layers in the vicinity of solid walls. When the Prandtl number 
is large we expect that thermal boundary layers will be much thinner 
than viscous boundary layers. And when Damkohler numbers are large we 
expect narrow reaction fronts. 

In chemical kinetics we find large numbers of chemical reactions com- 
peting for the available chemical species at widely different rates. As a 
result some species are either consumed or produced rapidly or slowly, 
while other species are being both produced and consumed simultaneously 
at high rates, with their net production rate being relatively slow. This 
chemical behavior is responsible for the many widely differing time scales 
in the mathematical models. The computational models of these processes 
are characterized as either multirate problems or stiff problems. 

It is worthwhile to point out the distinction between multirate prob- 
lems and stiff problems. In both problems the system itself is equally 
capable of rapid or slow changes. Multirate problems are those in which 
at least one component of the solution is changing rapidly, even though 
others are changing slowly. Numerical methods for these problems must 
take time steps that are small enough to resolve the fast transients, so they 
are controlled by accuracy not stability considerations. Stiff problems, on 
the other hand, are those in which all components of the solution axe 
changing slowly compared with the fastest characteristic scales possible in 
the model. In these cases explicit numerical methods are forced to take 
much smaller tune steps than are needed to maintain accuracy in order 
tomamtam stability. Often problems that begin as multirate problems 
become stiff problems as an equilibrium or steady-state condition is ap- 
proached. Stiff problems are usually solved efficiently by implicit methods. 
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For stiff or multirate problems, it is perhaps useful to consider methods 
that treat the fast components differently than the slow components. 

Given that multiple-scale problems are of practical importance, we 
must consider: why are there computational difficulties, and what can be 
done to ameliorate the difficulties? In the case of disparate length scales 
one difficulty is that of representing accurately the highly localized be- 
havior. If the solution is represented discretely on a mesh network, the 
mesh must be sufficiently fine to capture the localized behavior accurately. 
The whole topic of adaptive meshing is critically important for these prob- 
lems (see Section 4.2). Here, instead of computing on a fixed prespecified 
mesh, the mesh adjusts itself dynamically as the solution develops hi order 
to maintain accuracy in the solution. 

For situations in which the localized behavior is known to be ap- 
proximated well by very sharp fronts (e.g., shocks or flames), front track- 
ing methods can have significant advantages. Unlike the adaptive mesh 
approach where the solution is resolved smoothly through the front, the 
front tracking methods approximate the front by a discontinuity whose 
magnitude, speed, and location are to be found. Then elsewhere in the 
region the conventional discrete representations are adequate. 

In problems like chemical kinetics, the disparate time scales cause 
the governing differential equations to be stiff. Here, explicit solution 
methods are well known to be extremely inefficient, and some form of 
implicit method is needed. For systems of ordinary differential equations 
the problem has been worked out, and high-quality computer software is 
available. However, when stiffness is encountered in the context of systems 
of partial differential equations the remedies are much less developed. The 
same techniques used for ordinary differential equations, when applied 
directly the partial differential equation problems of practical interest, 
often yield problems that are simply too large for current computers. 

Several approaches show promise. One is to develop operator- splitting 
methods in which the stiff parts of the problem are split off and solved as 
a series of smaller and hence more tractable problems. Another approach 
is to attempt to remove the stiffness by solving instead an approximate 
(yet sufficiently accurate) system of equations. This tack benefits from an 
asymptotic analysis of the equations. Both approaches have found recent 
successes in fluid mechanics and in combustion chemistry. 
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3.3 SINGULARITIES IN COEFFICIENTS, DATA, OR STATES 

Difficulties similar to those encountered in multiple-scale models are often 
found in problems having singularities in coefficients or states. That is, a 
singular or discontinuous coefficient can give rise to localized behavior in 
the solution, such as very steep fronts. An example could be a material 
interface in a structural or heat-transfer problem, say between a steel and a 
plastic part. At this material interface the solution (stress or temperature 
gradient) might change rapidly. In order to maintain accuracy in the 
computed solution, the numerical procedure would have to resolve this 
frontlike behavior. The situation is analogous to the occurrence of a shock 
or a flame front. However, we usually know where the material interfaces 
are, so they are perhaps easier to deal with than phenomena such as shocks. 
Phase transitions can also produce discontinuous coefficients. Take 
a model in which a melting front is traveling through a region. Usually 
the properties of the molten material are quite different than those of the 
solid material. In fact, different governing equations may even be required 
for the two regions (e.g., fluid motion may be modeled in the liquid but 
not the solid). In any case, the solution is likely to exhibit jumps in its 
properties at the transition, and the numerical method will need to locate 
and resolve it. This situation is more like a shock, in that the position of 
the phase transition front is not known a priori, and thus the numerical 
method must both locate and resolve it. 



3.4 BOUNDARY CONDITIONS 

The solution to a boundary or initial-boundary-value problem depends 
strongly on the boundary conditions. Thus it is important to understand 
the relationship of the boundary conditions to the differential equations 
and to their discrete representation. Most important, the boundary con- 
ditions must be chosen so that the problem is well posed. For a large class 
of problems there is a satisfactory theory of admissible boundary condi- 
tions, but for many problems, those involving coupled hyperbolic-elliptic 
systems or disparate time scales, for example, only a rudimentary theory 
is available. 

A possible error in prescribing boundary conditions for hyperbolic 
equations is to overspecify or underspecify the number of boundary condi- 
tions. Overspecification usually causes nonsmooth solutions with mesh os- 
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dilations near the boundary. Underspecification does not ensure a unique 
solution, and the numerical solution may tend to wander in steady-state 
calculations. In either case, the results are not accurate and one should 
be skeptical of even the qualitative behavior of the solution. It should 
be noted that the way in which boundary conditions are specified for the 
difference equations can change a well-posed continuous problem into an 
ill-posed (unstable) discrete problem. 

Two of the most common methods used to incorporate boundary 
conditions into discrete equations are the extrapolation and uncentered 
differences methods. In the extrapolation method, the domain of the prob- 
lem is extended and the solution is extrapolated to fictitious points outside 
the integration region. The nonphysical solution at these points is defined 
so that the discrete equations are consistent with as many relationships 
as can be derived from the boundary conditions and differential equa- 
tions. The extrapolation formula can do this best by incorporating the 
discrete boundary conditions into the extrapolant. Additional relations 
can be generated by differentiating the boundary conditions with respect to 
time, replacing all time derivatives by space derivatives using the original 
differential equation, and discretizing the resulting equations. 

The uncentered differences approach is to extend the number of 
boundary conditions so that all components of the solution are defined 
at the boundary. Again, these additional boundary conditions must be 
consistent with the original problem and as many relationships as can be 
derived from it. An uncentered difference approximation is then used to 
approximate the spatial derivatives at the mesh points nearest the bound- 
ary. 

Irregular domains can be imbedded in an underlying regular grid 
that is not aligned with the boundary, or an irregular grid can be con- 
structed that conforms to the boundary. The discrete approximations to 
the equations away from the boundary are much simpler on the regular 
imbedded grids, but the boundary conditions are difficult to approximate. 
Boundary-fitted grids can be generated algebraically in the original physi- 
cal domain, or the domain (and hence the grid) can be mapped onto a 
regular grid in a simpler domain and the equations solved there. The 
algebraic-grid-generation methods have the advantage that the equations 
and boundary conditions are unchanged, but the differential operators are 
more difficult to approximate on the nonuniform grid. When using the 
mapping method, the differential operations are easily approximated but 
the transformation can greatly complicate the equation and sometimes 
obscure important properties such as the conservation laws expressed by 
the equations. 
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3.4.1 Boundary Conditions at Infinity 

Many physical problems require the solution of partial differential equar 
tions on some infinitely large domain fL For computational reasons this 
domain is often replaced by a finite domain QI. Then the difficult problem 
of specifying boundary conditions at its finite artificial boundary B arises. 
It is especially important that these artificial boundary conditions do not 
introduce spurious phenomena. Consider, for example, a nonviscous fluid 
that at subsonic speed leaves HI through the boundary B. There is one 
characteristic direction that points back into the region HI, and there- 
fore one boundary condition has to be specified on B. But, in general, 
no detailed knowledge of the solution on B is known and therefore other 
principles have to be applied. For example, if one has solved the problem 
by difference approximation then one predicts the solution on B from in- 
side wi, for all the dependent variables. Thus this procedure amounts to 
overspecification of the solution on B. Another principle has been proposed, 
namely, to specify the boundary conditions on B so that no reflection of 
high frequency takes place. However, numerical experiments have shown 
that such approaches do not always work. 



3.5 BIFURCATIONS AND CHAOS 



3.5.1 Buckling and Collapse Behavior, Bifurcations 

In general, the equilibrium equations of a mechanical structure involve a 
finite number of parameters, that is, they have the generic form F(x, p) = 
0, where x varies in some state space X and p 6 R m represents a parameter 
vector. Thus, in general, the solution set {(x,p) 6 X xR m ; F(x,p) = 0} is 
a manifold X x/?" 1 , and one topic of interest is the location and character 
of the singular points of the solution set. 

Without going into detail, suppose that we follow some curve on 
the solution manifold defined by some combination of parameter values 
with one degree of freedom represented by a single variable X. Then 
we may encounter certain critical points on the solution curve where 
the mechanical structure may suffer a loss of stability. Such a loss of 
stability actually corresponds to a dynamic phenomenon whereby the 
structure undergoes a sudden change of deformation. The dynamics of this 
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phenomenon are not described by the equations of equilibrium F(x, p) = 0, 
but it is possible to deduce from the shape of the (static) solution manifold 
at that point what type of sudden changes may be expected. 

We use a few figures to illustrate the situation. In the upper left-hand 
part of Fig. 3.1 the point denoted by (1) is a so-called limit point or turning 
point and an increase of the load-intensity X beyond the critical value X c 
results in a jump from (1) to (presumably) (2). This type of behavior is 
called a snap-through or collapse. In the case of the upper right-hand part 
of Fig. 3.1 the instability phenomenon is related to the bifurcation of the 
solutions at (1). This behavior is called buckling. This part of the figure 
is a classical example of stable buckling where a distinct change in the 
character of the load deformation is encountered when the load-intensity 
X passes the buckling load X c but where no failure of the structure occurs. 
On the other hand, the lower part of Fig. 3.1 shows an unstable buckling 
point where we again observe a dynamic departure from the bifurcation 
point to some other state [presumably the equilibrium state (2)]. The 
geometrical shape of the bifurcation branch n is the determining factor in 
the question of whether the bifurcation point is stable or not. 

These examples already indicate that for a deeper understanding of 
the behavior of a mechanical system it is necessary to analyze the shape of 
the full solution manifold. Of course, the choice of the parameters entering 
into the definition of the equation is of critical importance here. In essence, 
catastrophe theory provides some information about the selection of par- 
ticular minimal numbers of parameters, but in practice the parameters 
are simply chosen to correspond to the certain natural features of the 
structural problems. 

In view of these observations the aim is to develop procedures for 
a computational analysis of the form of the equilibrium surface. Some 
methods for this purpose are mentioned in Section 4.4 on continuation 
methods, but the entire problem is still a wide open research problem. In 
particular, it has to be noted that we can compute only points that are 
approximately on the solution manifold of a discretization of the original 
problem. Thus, the questions arise whether phenomena, such as limit 
points, or stable or unstable bifurcation points, encountered on the solution 
set of the discretized problem actually correspond to similar phenomena 
for the original problem, and, if so, what errors have been encountered. 
These questions are as yet largely unanswered and represent considerable 
research challenges. In this connection, it might be mentioned that for 
nonlinear problems of this type the solution manifold of the discretized 
equations often has a different number of connected components than 
that of the original equations. The components that do not approximate 
the exact solution manifold have been called spurious, or numerically 
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FIGURE 3.1 Instability behavior of shells under different loadings F. On 
the right the load intensity X is plotted versus a characteristic deflection 
W. [Prom W. Rheinboldt and E. Riks, "A survey of solid techniques 
for nonlinear finite element equations," State of the Art Surveys of Finite 
Element Technology, Chap. 16, special AMD pub!., American Society of 
Mechanical Engineers, New York (1983).] 
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irrelevant solutions. They are being observed more and more often in 
various applications, but their study has only recently begun. 



3.5.2 Chaos in Deterministic Systems 

It has been remarked that one of the most engaging problems in nonlinear 
dynamics is that of understanding how simple deterministic equations can 
yield apparently random solutions. As is now widely recognized, systems 
with this behavior appear hi many of the branches of scientific endeavor. 
This near universality has come to the attention of physical scientists and 
applied mathematicians principally through the use of computers in the 
study of properties of dynamical systems. 

The onset of chaos in deterministic systems (the stochastic instability) 
signals an unusual sensitivity to initial conditions, i.e., the trajectories hi 
phase space diverge widely as tune goes on even though the initial con- 
ditions are arbitrarily close to one another. The behavior is such as one 
would expect hi a space with everywhere negative Gaussian curvature. On 
the other hand the trajectories may hi fact tend to a single orbit nearly 
filling a subspace of the phase space (a strange attractor) in an ergodic 
manner. In even the simplest systems of this type it is possible for several 
such attractors to coexist side by side, with the initial conditions deter- 
mining which one is reached asymptotically hi tune. These and other 
numerous unusual properties of dynamical systems displaying chaotic be- 
havior have led to some open problems in computational mathematics. 

Examples of such problems include the following: When on varying a 
parameter the stochastic instability sets in, the continued use of an algo- 
rithm describing the evolution of the system prior to the onset of chaotic 
behavior may no longer be appropriate. In this event it may be more prac- 
tical to take advantage of the stochastic nature of the system and use more 
or less conventional statistical methods. In order to effect such a change in 
computational methods it is necessary to detect the change from orderly 
behavior of the system to chaotic behavior. Theoretically, the sensitivity 
of the system to initial conditions accompanies the occurrence of positive 
Lyapounov exponents, i.e., numerical indices of the asymptotic rate of 
divergence of initially nearby system trajectories. Since the Lyapounov 
exponents are defined as time asymptotic quantities, straightforward com- 
putation of these quantities requires a considerable computational effort, 
including a calculation of the evolution of two initially nearby trajectories, 
or the simultaneous integration of the associated variational equations. 

An important tool for studying the properties of dynamical systems 
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in the chaotic regime is the first return map, sometimes referred to as 
the Poincar6 map. Such a map displays in a graphical manner various im- 
portant aspects of the attractor. Unfortunately, the simplicity of generat- 
ing such a map is restricted to systems no larger than three dimensional, 
because of the self-evident difficulty of making multidimensional graphic 
displays. As the study of higher-dimensional dynamical systems advances 
there -will be an ever more urgent need for a higher-dimensional equivalent 
of the first return map. 

As mentioned earlier, in the chaotic regime it may be more appropriate 
to describe the properties of the system trajectory in statistical terms than 
in terms of a trajectory evolution. In order to obtain such a statistical 
description it is necessary to have appropriate information about the in- 
variant measure associated with the given dynamical system. With few 
exceptions such measures cannot be derived a priori but must be obtained 
from detailed calculations of the trajectories. A finite machine computa- 
tion of the system trajectory will inevitably introduce some errors that 
may be all the more serious in the chaotic regime because of the previously 
mentioned exponential divergence of nearby trajectories. There is then the 
question of how accurate the computation of a chaotic trajectory must 
be in order to yield enough information for constructing the appropriate 
invariant measure. 

While the above examples have been drawn from dynamical sys- 
tems representable by ordinary differential equations there are other sys- 
tems, e.g., those described by partial differential equations and by integral 
equations, whose study is likely to be replete with similar computational 
difficulties. 



3.5.3 Symmetry Breaking 

Bifurcations commonly arise in connection with a loss, or breaking of 
symmetry. In such cases the extra structure of the problem symmetry 
may simplify the analysis. We give an example from fluid mechanics, A 
classic experiment concerns flow between two rotating vertical cylinders. 
For small angular velocities of the cylinders the flow is laminar and is called 
Couette flow. As the angular velocity is increased the vertical translation 
symmetry is broken and axisymmetric Taylor cells appear. These cells 
resemble a row of adjacent smoke rings, each one rotating in a sense 
opposite to its immediate neighbor. A further increase in velocity breaks 
the axial symmetry, and the Taylor cells become wavy with a time periodic 
shape. Eventually the flow becomes turbulent. Specific experiments on 
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these systems are a classic topic in fluid mechanics. However as often 
happens in science, one cannot see important phenomena in the absence of 
a carefully thought out theory. Recently a number of precise experiments 
have been conducted concerning the role of strange attractors. Routes to 
turbulence and, hi particular, wavy Taylor cells were carefully observed. 
The results confirmed some of the theoretical predictions but not all of 
them. One theoretical proposal was that generic limiting sets, known as 
strange attractors, would dominate turbulent flow patterns. This proposal 
was supported by the fact that strange attractors were observed in low- 
mode approximations to fluid flows. More refined calculations with more 
modes included indicate an absence of these strange attractors, but the 
phenomenon is still indicative of the highly complex solution manifolds 
that can arise hi nonlinear problems. 



3.6 ILL-POSED INVERSE PROBLEMS 

The notion of a well-posed problem is due to Hadamard: a solution must 
exist, be unique, and depend continuously on the data. The term "data" 
can have a variety of meanings; in a differential equation it could include 
any or all of the following: boundary values, initial values, forcing term, 
and even the coefficients in the differential equation. Since data cannot be 
known or measured with arbitrary precision, it was felt for a long time that 
real physical phenomena had to be modeled by mathematically well-posed 
problems. This attitude has changed considerably in recent years, and it 
is now recognized that many applied problems are ill-posed, particularly 
when they require numerical answers in the presence of contamination of 
the data. 

El-posed problems often arise in the inversion of well-posed problems. 
Consider, for instance, a well-posed problem that smooths the data or 
attenuates its high frequencies. The inverse problem, in which the role of 
data and solution are interchanged, will then be ill-posed. A simple but 
important example is the Fredholm integral equation of the first kind 



/ k(x, y)u(y)dy = v(x), < x < 1 
Jo 

or, in operator form 

Ku = v 

Assuming the kernel to be continuous on the closed unit square, the 
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Riemann-Lebesgue lemma gives 

yi 
lim / k(x t 2/)sin nydy = 

n ^oo JO 

so that K attenuates high frequencies and thus transforms widely different 
functions u into approximately the same v. If v is regarded as the data, the 
solution u of the integral equation, therefore, cannot depend continuously 
on the data. Indeed, it is not necessarily uniquely solvable, either take, 
for example, 
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If, unaware of these difficulties, one attempted to solve the integral equa- 
tion by discretization, one would find that the corresponding matrix prob- 
lem is singular or ill-conditioned (and the finer the discretization, the more 
ill-conditioned the matrix problem). 

More general inversion problems can also be reduced to an equation 
of the form Ku = v, where K is a continuous transformation, v is the 
data, and u is the solution being sought. It may happen that u itself is 
not the quantity of principal interest but rather some functionals of u such 
as some of its moments or its values at a few specified points. 

Thus, there are three pieces of information that are central to the 
numerical resolution of an inverse problem: (1) the model M, representing 
the equation involving a mapping K between appropriate spaces; (2) the 
observation operator 0, representing the measurements that can be made; 
for instance, we might have O(w) = {v(xi)-,i = l,2,.,.,n}; (3) the intel- 
ligence operator J, which specifies the information we wish to extract from 
the solution. For instance, we might have J(u) = {fx k u(x)dx; k = 0, 1, 2}. 

Most of the existing approaches focus almost exclusively on the model 
M, taking account of a priori information about the solution such as 
smoothness, positivity, and bandwidth. The Tychonov regularization 
method is of this kind. By restricting K to a suitable subspace it may 
be possible for the restriction of AT to be one to one with a continuous in- 
verse. It might be more appropriate to study the triplet M,0, J for different 
choices of O and J. One goal of such a new approach would be to classify 
and quantify ill-posedness by developing comparison principles and order 
structure. 

Hi-posed problems with their attendant numerical difficulties abound 
in practice. For instance, in scattering from obstacles one may wish to 
determine the shape of an object (or its surface impedance, if its shape is 
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known) from far-field measurements. Coefficient identification problems 
arise in many contexts, one such problem being the determination of the 
sound speed in subsurface media from measurements of the field at or 
near the surface. Another class of ill-posed problems is associated with 
image reconstruction (tomography), the medical applications of which are 
now widely used (CAT scanners, for example) (see Section 2.12). More 
recently, similar ideas have been applied to the nondestructive testing of 
mechanical structures to detect cracks in fuel rods, weaknesses in rotor 
blades of jet engines, and faults in screws in ships. Optimum filtering 
and inverse problems in Fourier optics (restoring data designed by a band- 
limited filter, for example) are other areas of current research. 



3.7 EFFECTIVE MEDIA 

The study of "bulk" or "effective" parameters for composite media is of 
fundamental importance. Depending on the particular application area it 
may be relevant to consider a periodic (or an almost periodic) or a random 
formulation. 



3.7.1 Homogenization, the Deterministic Approach 

Many problems of physical interest involve several length scales. As an 
important example we mention the study of composite materials in struc- 
tural mechanics. Owing to a particular manufacturing process, a distinct 
structure, e.g., periodically or almost periodicity, is often present. 

Homogenization is an approach for deriving the macroscopic 
properties of the material from the known microscopic ones. A variation 
of this is to replace a complicated geometrical configuration by a simpler 
one, e.g., to replace a framework by a plate or to smooth out a rough 
surface. 

Homogenization may be applied to linear and nonlinear problems and 
can provide qualitative information about physical macrolaws as well as 
good approximations to the various parameters appearing in these laws. 

The most important mathematical tool for the deterministic approach 
is some form of asymptotic expansion. The theoretical results concern the 
limiting behavior when certain of the relevant scales become very small. 
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One is likely to obtain different homogenized formulations from different 
limiting relations between the length scales. 

The numerical solution of a problem introduces a new length scale: 
the mesh size in a finite-element method, the step size in a finite-difference 
method, or the wavelength in a spectral method. Different limiting rela- 
tions between this new length scale and the original physical length scales 
usually lead to different results. In this sense the numerical treatment has 
to be considered simultaneously with the homogenization process. One has 
to design algorithms that will adaptively select the correct homogenized 
formulation and discretize it appropriately. 

Whereas much theoretical work has been done within the last ten 
years analyzing the effects of multiple scales in continuous media, the 
study of numerical discretizations of such problems is still very much in 
its infancy. 



3.7.2 Analysis of Random Media 

Assessing the effect of random fluctuations in the coefficients of a par- 
tial differential equation is a basic mathematical problem that arises con- 
stantly in science, applied science, and engineering. Finding effective con- 
ductivities of composite conducing materials such as soil or a metal alloy, 
finding effective fluid equations for flow in porous media (Darcy's law), 
determining the rate of sedimentation of particles in a fluid, and many 
other problems are problems that may require analysis of random media. 
In addition to the quantitative aspects of the problem, many interesting 
qualitative questions can be posed as well. For example, what is the nature 
of the spectrum of the Schrodinger equation with a random potential? It 
has been shown that in one dimension, for a large class of random poten- 
tials, the spectrum is always discrete. For randomness of large magnitude, 
it was recently shown that there is no diffusion, without restriction on the 
dimension of the space. 

From the applications' viewpoint, one can frequently model ade- 
quately a random medium, such as a suspension, for example by a con- 
tinuum with suitable constitutive properties. The continuum equations 
may be linear or nonlinear, and the constitutive laws may be known 
only qualitatively from experimental data. In such a context, theoreti- 
cal investigations are useful and necessary in order to understand how 
the phenomenological continuum equations arise from the known micro- 
scopic Jfrueta* One can then find mathematical characterizations for 
the relevant constitutive laws that can lead to interesting conclusions. 
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Examples of this arise frequently as, for instance, the determination of 
bounds for the effective conductivity of composites. 

The effective media ideas of Maxwell in the last century dominated 
theoretical calculations for a long time. In the last two decades, press- 
ing technological needs have caused a major expansion in materials re- 
search. In attempting to achieve a mathematical understanding of the 
methodological basis for effective media calculations that abound, one finds 
a lack of theoretical foundations. 

Mathematical methods in random media are drawn from analysis, 
probability theory, asymptotic methods, and differential equations. The 
goal is to develop tools that bridge the gap between microscopic and 
macroscopic descriptions, give qualitative information about constitutive 
laws, and determine when residual stochastic effects remain and how they 
can be characterized. In addition, it is of interest to find important specific 
problems on which more detailed analysis, including numerical analysis, 
can be carried out so that these problems can serve as benchmarks in the 
theoretical development. 



3.8 VALIDATION, ERROR ASSESSMENTS, AND 
SENSITIVITY ANALYSIS 

The results look good, but are they really right? The question is often 
not answered satisfactorily. Indeed, often too little attention is paid to 
the difficult topic of model and code validation. Once the results look 
plausible, we are often either unable or unprepared to take the valida- 
tion process further. The validity of a model depends both on having a 
proper physical model (Do the governing equations adaquately represent 
the physics?) and on having an accurate computational representation of 
the governing equations. The mathematician has a responsibility in both 
areas. He should help determine the "well-posedness" of the models and 
from a mathematical point of view help the physical scientist determine 
the appropriateness of the model. Given that the equations are proper, the 
computational mathematician must be sure that the numerical procedures 
used accurately approximate the solutions of the governing equations. 

Frequently modeling is done in conjunction with experiments, and 
those results are used to validate the model. Comparison with experiments 
simultaneously tests both the validity of the governing equations in the 
model and of their numerical solution. In this case the source of any 
discrepancy is not easily isolated. A reasonable goal should be to validate 
the numerical procedures and the physical models separately so that a 
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model can be developed more confidently. The design of reliable numerical 
error estimates for the computational methods is an important step. 



3.8.1 A Posteriori Error Estimates 

An important current research topic is the development and application 
of a posteriori error estimates. Such procedures are valuable for models in 
which differential equations are approximated by their discrete analogs on 
a mesh network. When a posteriori estimates of the error associated with 
the discretization of the continuous model are possible and they often 
are then potentially it is possible to control the errors. Adaptively moving 
the mesh network to control or mimimize the error is one application of 
this approach. With a posteriori error estimates, it should be possible, 
in principle, to give a strict bound for the error on completion of the 
computation. 

The selection of the specific accuracy requirements depends strongly 
on the goal of the computation. Often it is desired to obtain detailed 
information about the solution itself; in other cases, the main focus is the 
value of a specified functional of the solution, as, for example, a stress in- 
tensity factor in fracture mechanics or a drag coefficient in fluid dynamics. 
Other goals may be the determination of certain critical data, such as 
collapse points or buckling points and their associated loads in structural 
mechanics. Sometimes special techniques, e.g., variational techniques, can 
yield highly accurate estimates of particular functionals. 

In connection with most of these goals we are interested in quantita- 
tive results that have a desired accuracy. For this the error has to be 
defined, in that a family of exact results has to be specified with which 
the computed data are to be compared, a norm has to be prescribed in 
terms of which the error is to be measured, and some procedure has to be 
established for estimating the error. 

Such error estimation capabilities are certainly important hi many 
applications, provided they can be guaranteed to be reliable. For example, 
for the many types of certification computations required hi the design of 
complex structures or nuclear plants, the availability of reliable estimates 
of the accuracy of the computed data is obviously important. In other 
cases, such estimates may reduce the total design effort and avoid un- 
necessary overdesign. At the same time, the availability of effective er- 
ror estimates introduces the possibility of applying adaptive techniques to 
structure the computation to achieve a desired error tolerance at minimal 
cost or to provide the best possible solution within an allowable cost range. 

Many of the theoretical studies of solution algorithms for classes of 
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differential equations provide for some a priori error estimates, usually 
of an asymptotic type. For the practical error assessment these a priori 
bounds are rarely computable or very accurate. Thus, one is led to the 
necessity of developing a posteriori error estimates that utilize information 
gained during the computation itself. 

For the finite-element solution of certain classes of elliptic boundary- 
value problems, some computable and reliable a posteriori error estimates 
have been developed and analyzed in recent years. Most of these a pos- 
teriori estimates are based on a local analysis. For a given mesh A consist- 
ing of elements AI, . . ., A n and with a corresponding approximate solution 
on that mesh, an error indicator rjj is associated with the yth element A J( 
These r\j have to be computable in terms of information about the prob- 
lem and the approximate solution on Aj and, at most, on the immediate 
neighbors of that element. On the basis of the indicators r}i,...,r) n an 
error estimate e(A) is then constructed. Of course, the T)J as well as e(A) 
depend on the chosen norm. The theoretically important question is then 
the relationship between e(A) and the norm of the actual error e(A). The 
effectivity of the estimation may be judged by the effectivity index 9 = 
(A)/||e(A)||. In practice, it is usually more important for 6 to be close 
to 1 than that < 1. Moreover, it is essential that 6 converges to 1 when 
||e(A)|| tends to zero, so that for an accuracy of, say, 5 to 10 percent the 
value of \0 1| is expected to be less than, say, 0.1 or 0.2. 

The requirement of designing error estimators with these properties 
for realistic classes of problems and various different norms certainly repre- 
sents a demanding research task. The results available so far suggest that 
such error estimators can be developed at least for linear problems. For 
nonlinear problems the situation is in an embryonic stage although some 
results available for model problems in one space dimension indicate that 
estimators for the error along continuation paths and for the location of 
critical points are computationally feasible. There is certainly considerable 
need for concentrated research in the general area, 



3.8.2 Other Validation Measures 

Checking for mesh dependencies and convergence rates ie another way to 
help validate a model. If a method is supposed to be, say, "second order," 
then as the mesh intervals are halved, the errors should be reduced by a 
factor of 4. If this does not happen according to one's error estimate, then 
the method does not have the desired order yet and may be in error. Also 
one should continue to refine the mesh until the results are independent of 
refinement to within the desired accuracy. Solutions can also be checked 
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for sensitivity to mesh orientation. The solution should not be aligning 
itself with some special property of the mesh network. 

Numerical methods can also be validated by comparing them with 
analytic solutions. Often this is possible by choosing limiting cases where 
an asymptotic analysis applies. In some cases complex constitutive laws 
can be replaced with constants so an exact solution is possible. Higher- 
dimensional models can sometimes be validated by comparison with pre- 
viously confirmed lower-dimensional results. Bounds on solutions are 
another way to help verify a model. Often conservation laws or minimum 
and maximum principles apply, or the solution is known to approach a 
previously known equilibrium or steady-state condition. These properties 
can be used as a basis for validation. 

To the extent possible one should always write computer codes in 
a modular fashion such that each part of the model can be validated 
separately. Even though this is just good structured programming prac- 
tice, and not necessarily related to mathematics, it is an important aspect 
of model validation. 



3.8.3 Sensitivity Analysis 

In addition to the independent and dependent variables, most models also 
depend on certain physical parameters. Unfortunately, those parameters 
are often not known accurately. Sensitivity analysis is a systematic 
means to quantify the relationship of model parameters and model results. 
Doing the sensitivity analysis requires solution of an additional system 
of differential equations. These equations are a formal statement of 
the relationship between the dependent variables of the system and the 
parameters. The results are given in terms of a matrix of partial deriva- 
tives of the dependent variables with respect to the system parameters. 
One way to think of the analysis is as a method to provide theoretical 
"error bars" for the model. 

New methods are being developed to solve these sensitivity equations 
quite efficiently. The methods rely on the observation that the sensitivity 
equations are linear, regardless of the nonlinearities in the original model. 
The method solves differential equations for the Green's function of the 
sensitivity equations, and then the sensitivity coefficient matrices are com- 
puted by quadrature for the various imhomogenous terms corresponding 
to the parameters. These methods have been successfully applied to prob- 
lems in chemical kinetics. 
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4. NUMERICAL METHODS 



4.1 DISCRETIZATION METHODS 

Many important physical problems are modeled by boundary or initial- 
boundary-value problems. In this modeling, the physical state under 
consideration is characterized by a function u, which is the unique solution 
to the boundary or initial-boundary-value problem. Thus a major part of 
the analysis of the physical problem is the determination of u or some of 
its derivatives or some functional of it. 

In nearly all important problems the determination of u is an infinite- 
dimensional problem in the sense that u does not lie in an explicitly 
known finite-dimensional space or, alternatively, cannot be expressed in 
terms of a set of explicitly known functions by means of a finite set of 
parameters. Thus u, cannot be "completely" represented on a computer, 
and it is necessary to resort to some type of approximation or discretization 
method. In essentially all discretization methods one attempts to construct 
a function, u &pprox , which on the one hand is "close" to u and on the other 
is characterized by a finite set of parameters, and so can be calculated. 

The aim of finding ii approx that is "close" to u is made precise by 
requiring that \\u - u &pprox \\ < r (or that \\u - u appr ox || < r||u||), where 
1 1|| is a physically relevant norm and r is a physically relevant tolerance. 
The goal is thus to find u app rox satisfying this error criterion with the 
least expenditure of computational effort. The selection of a discretization 
procedure is influenced by a number of considerations, the most important 
of which we now list: 

1. The goals of the analysis of the physical problem. 

2. Known mathematical properties of the physical problem and the 
algorithm. 

3. Hardware considerations, e.g., the availability of parallel process- 
ing. 

4. Computer-science considerations, e.g., data-management require- 
ments, 

5. Restrictions on computation time and expense. 
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As indicated above, the exact solution is approximated by a function 
u a pprox , which is expressible in terms of a set of explicitly known (basis) 
functions by means of a finite set of parameters, which are determined 
in the computation. One broad classification of methods is in terms 
of the nature of the basis functions, namely (a) those involving global 
basis functions having global support and (b) those involving local basis 
functions having small support. Another classification is based on tne 
extent to which the method is adaptive. Adaptive methods refine and 
modify themselves on the basis of partially completed computations. 

Of course, these classifications are not precise, and there are methods 
possessing the different properties in various degrees. We now turn to a 
more detailed discussion of the most commonly used discretization proce- 
dures. 



4.1.1 Finite Differences 

One of the most frequently used discretization methods is the method 
of finite differences. The central idea is to replace each partial deriva- 
tive occurring in the differential equation in the underlying boundary or 
initial-boundary-value problem by an approximating difference quotient, 
For example, the first-order derivative du(x,y)/dx may be replaced by 
the forward difference [u(x + h,y) -u(x,y)]/h or the backward difference 
[u(x,y] - u(x - fr,y)]//i, whereas the Laplacian Au = u xx + Uyy may be 
replaced by the five-point difference operator 



u(x + h.y} + u(x. v + h) + u(x- h. 



This replacement of all derivatives in the differential equation by ap- 
propriate difference quotients leads to a system of equations, called 
difference equations, for the numbers Uij, which are to approximate the 
values of the exact solution u(x,y] at the finite-difference mesh points 
(yhi kh), j, k = 0, 1, +2, . . ., where h is a small positive number, called the 
mesh parameter. Nonuniform meshes may be considered as well, but at 
considerable complication in the form of the difference quotients in the 
resulting difference equations. 

There are a number of important questions that arise in connection 
with the study of difference methods, namely questions of 
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1. Accuracy of the method; 

2. Stability of the method; 

3. Efficient solvability of the finite-difference equations; 

4. Robustness of the method with respect to the input data (e.g., 
coefficients, forcing functions, meshes); and 

5. For nonlinear problems, the separation of legitimate approximate 
solutions, which correspond to an exact solution, from spurious ap- 
proximate solutions, which do not correspond to any exact solution. 



Finite-difference methods are an example of the local basis class of 
approximation methods. They can be adaptive. Adaptivity is usually 
introduced via adaptive mesh selection: the mesh chosen at any stage in 
the computation is based on the previous computations. 

Much important work on difference methods remains to be done. We 
mention in particular the construction of effective difference methods for 
nonlinear problems, in particular for those problems whose solutions have 
shocks. 



4.1.2 Variational Methods of Discretization 

We will now discuss a class of discretization methods based on variational 
formulations of the physical problem under consideration. As noted above, 
the problem of calculating the exact solution u is infinite dimensional in 
the sense that u is only known a priori to lie in an infinite dimensional 
space, say H. For a linear problem, a variational method of discretization 
consists of (a) a finite-dimensional space S cH, called the trial space, 
in which the approximate solution is sought, (b) a finite dimensional test 
space V, and (c) a bilinear form B(u, v) defined on H X V. The approximate 
solution is then determined by requiring that 



u 



approx 



-B(ii a pprox > i>) = B(u, 11), for all v 6 V 

where for veV , B(u, v) is computable from the data of the problem without 
knowing u, Since and V are finite dimensional, u app rox can De calculated 
by means of the solution of a system of linear equations if the original 
equation is linear. Usually, for nonlinear problems, B is nonlinear in u, 
and the resulting system of equations becomes nonlinear. 
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The approximate solution depends on the selection of S,V, and B. 
For any problem there exists a wide variety of variational methods of 
discretization, i.e., a wide variety of choices for S, V, and B. The rational 
selection of S, V, and B is a central problem. 

Variational methods can be of either the local-basis or global-basis 
type. Of those of global-basis type we mention the various versions of the 
spectral method and the p version of the finite-element method. Finite- 
element methods, considered in more detail below, and collocation methods 
are typically of the local-basis type. We also note that variational methods 
can be adaptive. 



4.1.3 Finite-Element Methods 

Finite-element methods arise if for 3 and V we choose spaces S h and V h 
of piecewise polynomials of fixed degree. That is to say, the underlying 
spatiaMomain for the problem is broken up into small geometric pieces, 
called "elements," whose size is measured by a parameter h. The functions 
in both S h and V h are then restricted to be polynomials in x and y on 
each piece but allowed to be different polynomials on the different pieces. 
In designing finite-element methods, i.e., in selecting S ht V h , and B, one 
attempts to achieve approzimability, stability, and systems of equations 
that can be solved effectively. 

Approximability here refers to the ability of the space S h to ap- 
proximate the unknown solution u. The solution u is unknown a priori, 
and often only the information u e H is available. In such situations 
6 h has to be selected so that every function in H can be approximated 
sufficiently well by one in S h . However, an approximation based on a few 
large elements can provide additional information on u, which can be used 
in turn to refine those elements. This type of adaptive element selection 
is especiaUy important for problems with sharply varying solutions (see 
section 4.2). 

In the choice of the bilinear form B one has, in effect, the freedom to 
choose a variety of variational principles, many of which have a natural 
physical connection with the original problem. For instance, the stan- 
dard Rita method is based on the principle of minimum potential energy. 
An alternate variational principle is the complementary energy principle. 
Approbation methods based on this principle are referred to as com- 
plementary energy or equilibribum methods. These methods involve a 
constraint that is usually difficult to satisfy. One way to circumvent this 
difficulty is to treat the constraint by means of a Lagrange multiplier. This 
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leads to the so-called mixed methods. They appear to be promising for 
many important problems and have recently received a large amount of 
attention. 



4.1.4 Transformation Methods 

In the transformation or pseudospectral method, the discrete approxima- 
tion u is first mapped by a transformation of the form 



into the m-dimensional function space of the coefficients, a,j. The basis 
functions fa and the transformation are chosen so that T and T~ l are 
fast (order of m log m operations) and so that differentiation D is simple 
in transform space. The derivative approximation can then be written as 



dx 

Some common transforms are based on the fast Fourier transform, where 
the fa are trigonometric functions or Chebyshev or Legendre polynomials. 
Selecting the <t> 3 - as piecewise polynomials with compact support, such as 
the B-splines, is another good choice. By choosing the transformation 
to incorporate some crucial property such as the periodicity or symmetry 
of u one can improve the accuracy of the method for a fixed number of 
basis functions. This can sometimes best be done by choosing the basis 
functions close to the eigenfunctions of the differential equation. 



4.1.5 Semidiscrete Methods 

When solving a partial differential equation one sometimes discretizes with 
respect to some but not all of the variables. For example, for the diffusion 
equation governing the cooling of a hot rod, one may discretize with respect 
to the space but not the tune variable, thereby replacing the original 
partial differential equation by a system of ordinary differential equations. 
Such an approximation method is referred to as a semidiscrete method or 
as the method of lines. Semidiscrete methods may be used for hyperbolic 
as well as parabolic equations. 
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One has, of course, eventually to discretize with respect to the time 
variable as well, Semidiscrete methods are based essentially on the view 
that very effective time discretization methods are available (in the form of 
preprogrammed software packages for the solution of ordinary differential 
equations) and that the spatial discretization is the main concern. An 
alternative point of view is to consider both discretizations simultaneously. 
Such fully discrete methods have been analyzed and tested extensively. 

Among important research problems for semidiscrete methods we 
mention the problem of adaptive mesh selection for the spatial discretiza- 
tion. 



4.1.6 Method of Characteristics 

This is a method for hyperbolic equations, particularly for those involving 
only one space dimension. In these equations the solution at some point 
in space-time depends primarily on its values along a fixed, finite number 
of curves (characteristics) going back in time from the given point. The 
approximations to these values are determined from difference equations 
that are closely related to the characteristics of the differential equations. 
This method has a natural generalization to quasi-linear second-order 
equations in two independent variables. It is especially important for 
problems whose solutions have shocks. 



4.2 ADAPTIVE GRID METHODS 

For realistic problems it is rarely feasible to design numerical processes that 
are rehable and accurate at a reasonable cost and yet that do not utilize 
some form of adaptivity. Put simply, most two-dimensional and virtually 
aJI three-dimensional problems are undercomputed without this technique. 
This statement will almost certainly remain true after the next generation 
of computers is available. The adaptive approach is to distribute the com- 
putational effort nonuniformly, so as to concentrate on the most singular 
parts of the solution or the most important parts of the solution from the 
point of view of critical design parameters. Correspondingly one must 
give less computational effort to the regular or less critical parts of the 

SL T 'I** SeCti ^ 4 ' 9) - M ** 8ame time > >-P**VB approaches 
may also lead to a simplification of the input data needed for the program 
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and hence free the user of part of the drudgery typical in preparing such 
input and in having to make the many a priori decisions required by most 
of today's programs. 

The goal of adaptive grid methods is to choose a grid that is par- 
ticularly refined, or that is aligned or oriented optimally, with respect to 
the solution in regions of space and time that are critical for solution ac- 
curacy. Thus adaptive grids utilize local mesh refinement or optimal local 
mesh orientation. 

The simplest adaptive grid is a preliminary, analytically determined 
coordinate transformation. For example, in the transformed coordinates 
the problem may be approximately independent of one variable or other- 
wise simplified. The next strategy is to choose a numerically deter- 
mined coordinate transformation. Typically, in two dimensions the coor- 
dinate transformation is obtained by the solution of a subsidiary partial 
differential equation. The resulting grids may be expected to give both 
improved mesh refinement and mesh orientation. The method is some- 
what problem dependent and can give rise to discretization errors in the 
mapping of solution variables between the various grids. 

Lagrangian grids for tune-dependent problems are adaptive for 
material interface singularities, because the grid points move with the 
material particles. Since these well-established methods also develop 
rezoning and mesh entanglement problems, they provide a reservoir of ex- 
perience concerning the difficulties associated with other evolving, adaptive 
algorithms. 

A refined grid can be constructed without recourse to a coordinate 
transformation. In response to a critical solution parameter or solution 
error criteria, selected regions of space can be flagged, preferred orientar 
tions selected, and refined subgrids introduced locally. Then, small time 
steps are chosen on the finer subgrid, and an interpolation problem must 
be solved to blend fine and coarse grid solution values. Finally, the con- 
struction is recursive, so that three, four, etc. levels of refinement can be 
introduced automatically, in response to some a posteriori error estimate 
on the next coarser grid. Precisely defined, reliable error estimators appear 
to be central to the design of effective adaptive procedures. 

In the context of finite elements for elliptic problems there is a devel- 
oped theory for a posteriori error estimates (and, hence, for adaptivity), 
which is based on local analysis. For a given mesh made up of elements 
AI, . . ., A n and with the aid of the corresponding approximate solution, an 
error indicator r)j is computed for each element Aj. For certain classes of 
problems it has been proven that asymptotically the errors become optimal 
when all the r\j become equal. This equilibration principle provides the 
basis for the control of the adaptive process. In essence only those elements 
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A.,- are refined for which the error indicators are too large in comparison 
with those of the other elements. And elements that are unnecessarily 
small are combined with their neighbors. The study of suitable algorithms 
for this has only recently begun. For example, some results have been ob- 
tained for algorithms based on the assumption that during the refinement 
process none of the elements could be combined without increasing the 
maximum value of the error indicators. 

In tune-dependent problems, differential equations can be introduced 
for the ever-changing optimal location of the grid nodes. These equations 
are then solved simultaneously with the original partial differential equa- 
tion, leading to what are known as moving mesh methods. 

Alignment of the grid without refinement is also possible, if fixed 
or moving boundaries or interior interfaces are specified as part of the 
problem. By alignment (without refinement), a regular grid index structure 
can be preserved. Maintaining the grid structure has the advantage that it 
potentially allows fast iterative methods, such as fast Fourier transforms, 
to be used as part of the solution algorithm. 

Adaptive grid techniques have been applied successfully to a range of 
time-independent problems. The newer time-dependent methods need to 
be developed to the point where they can be applied to meaningful two- 
and especially three-dimensional problems and compared with alternative 
methods. An important question requiring further attention is the con- 
struction of efficient, computable a posteriori error estimates for realistic 
classes of problems. Specifically, even for steady-state problems, there is 
the question of the validity of the equilibration principle mentioned above 
as well as the design and analysis of suitable adaptive control laws to 
implement this principle. This latter problem may require examination 
and incorporation of results in such fields as automatic control theory, 
artificial intelligence, and learning processes. It also raises the problem 
of the choice of appropriate refinement techniques that produce meshes 
with desirable properties. Last but not least there is the question of data 
management, which must be handled successfully to control the vastly in- 
creased algorithmic complexity associated with adaptivity and to achieve 
computational efficiency. 
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4.3 METHODS FOR SOLVING DISCRETIZED DIFFERENTIAL 
EQUATIONS 

Finite-difference and finite-element discretizations of partial differential 
equations usually give rise to large systems of equations in which each 
unknown is coupled to only a few of the other unknowns. Systems with 
tens or hundreds of thousands of equations are relatively common. For 
time-dependent problems these systems arise from the use of implicit 
time discretizations. For sufficiently fine grids, the numerical solution of 
these systems consumes a major part of the computer time for an entire 
simulation. 

Most nonlinear systems are solved by some form of iterative method 
based on linearization, such as Newton's method. At each step, these 
methods result in large sparse linear systems. In many cases, iterative 
methods converge only if the initial guess is sufficiently close to the solu- 
tion. 

There are two basic approaches to solving large sparse linear systems 
of equations: direct sparse matrix methods (i.e., some form of Gaussian 
elimination that takes advantage of the location of most of the zeros) 
and iterative methods, where the sparsity makes each iteration relatively 
inexpensive. No method is best for all problems. For many problems 
a combination of methods is attractive. Usually this takes the form of a 
block iterative scheme using a direct method to solve the subsystems whose 
diagonal blocks are matrices. Combination methods are of particular 
importance for the solution of linear systems arising from coupled systems 
of partial differential equations. 

Direct methods are relatively well understood today, and a number 
of excellent implementations of them exist for serial computers. Direct 
methods vary in the extent to which they take zeros into account. The 
simplest, nontrivial approach is band elimination, which takes account 
only of those zeros occurring outside the band determined by the extreme 
nonzeros of the matrix. The most complex approach is general sparse 
elimination in which all the entries that remain zero during the alimrnflr 
tion are taken into account. For systems that do not require pivoting 
(e.g., those having symmetric, positive definite matrices) there is a sym- 
bolic preprocessing phase in which the location of these zero entries is 
calculated. There exist good techniques (the nested dissection ordering 
and the minimum degree ordering) for arranging the unknowns and the 
equations so as to minimize the zero fill-in during elimination. For model 
problems it has been shown that the nested dissection ordering provides 
asymptotically optimal results with respect to work and storage. The 
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minimum degree algorithm is a valuable heuristic approach that is com- 
petitive in practice but that has defied analysis. For systems that require 
pivoting for numerical stability one cannot compute the zero structure a 
priori. Moreover, the ordering of the unknowns and equations to minimize 
zero fill-in mil usually be significantly altered. 

Some of the strong points of direct methods are as follows: (1) They 
provide an exact answer (modulo round-off error) to the linear system with 
a fixed number of operations independent of the system's condition num- 
ber. Most production structural analysis packages such as NASTRAN use 
some form of direct method, even for three-dimensional problems. In these 
applications (many of which are based on fourth-order elliptic problems) 
it is necessary to use higher-order precision because of conditioning prob- 
lems. (2) For problems with two dependent spatial variables, their execu- 
tion time is often shorter than the execution tune for iterative methods, 
especially for problems of moderate size. Some of their principal disad- 
vantages are the following: (1) they require considerably more storage than 
iterative methods, even for two dimensions; (2) they will almost always be 
noncompetitive with iterative methods for three-dimensional problems in 
terms of running time; and (3) except for the simpler methods such as 
band elimination (and to a lesser extent profile or envelope elimination) 
they do not vectorize well. In fact, owing to the necessary indirect ad- 
dressing involved, there are as yet no efficient implementations of general 
sparse direct matrix methods for the CRAY-1 or CDC CYBER-205 super- 
computers. 

Except for structural problems, iterative methods are commonly used. 
Classical iterative methods, such as successive overtaxation, Chebyshev 
semi-iterative methods, and such newer methods as the preconditioned 
conjugate gradient method, are fairly well understood for symmetric, posi- 
tive definite systems, and they are easy to implement. However, the situa- 
tion is not so bright for nonsymmetric or indefinite systems, though in 
practice classical iterative methods may converge rapidly with a clever 
choice of the iteration parameters. Nonsymmetric systems with indefinite 
symmetric parts are especially difficult to solve, and no satisfactory general 
algorithms are known at this time. Such problems arise in the simulation 
of secondary and tertiary thermal recovery techniques for petroleum reser- 
voirs (see Section 2.6). 

Some of the principal advantages of iterative methods are the follow- 
ing: 



1. They tend to require minimal storage (proportional to the number 
of unknowns); 
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2. They are reasonably fast for a wide range of problems. Moreover, 
the number of iterations required is independent of the number of space 
dimensions in the underlying partial differential equation (but not of their 
order, or of the mesh size); 

3. They can take advantage of good initial guesses to the solution, as 
would be available in time-dependent or nonlinear problems; and 

4. Many of them vectorize reasonably well on supercomputers. 



Some of their disadvantages follow: 

1. Mathematically rigorous stopping criteria may be difficult to for- 
mulate, e.g., for linear systems with matrices that are not symmetric, 
positive definite; 

2. Many of the methods require a selection of iteration parameters, 
and the performance of the methods depends critically on such a selection 
(this difficulty is being overcome somewhat by the relatively new adaptive 
methods and the preconditioned conjugate gradient-type methods); 

3. The interaction between linearization and iteration is not well 
understood, especially for discretizations of coupled systems of partial 
differential equations; and 

4. Nonsymmetric or indefinite problems may cause great difficulties 
for iterative methods. 



The relatively new multigrid iterative method combines the well- 
understood behavior of a given iterative technique with the fact that an 
underlying differential equation is being solved approximately. It alter- 
nates iterating, toward the differential equation's solution, on fine and 
coarser subgrids of the spatial network, with the goal of performing no 
more computational work on the finer (hence expensive) grids than is ab- 
solutely necessary. In many cases of practical interest, such as in neutron 
diffusion in complex environments, the technique yields sufficiently ac- 
curate solutions to the equations in a computational time proportional to 
the number of unknowns. This has been a long-sought-for goal in the 
approximation of elliptic equations in more than one space dimension. As 
an iterative method it also has a natural extension to nonlinear equations; 
and its logical structure, together with the already necessary calculation 
of residual errors, points toward the incorporation of adaptivity for the 
nesting spatial grids. 
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4.4 CONTINUATION AND HOMOTOPY METHODS 

For more than a century the so-called continuation principle has proved 
to be an important tool in mathematical theory. For example, early uses 
date back at least 100 years to H. A. Schwarz and his work on con-formal 
mappings; then, in the early part of this century, J. Hadamard and M 
Levy applied these techniques in connection with the inversion of nonlinear 
mappings, and it is also a basic tool in J. Leray and J. Schauder's work on 
differential equations. But it was essentially only in the early 1950s, with 
the advent of automatic computing, that continuation approaches began 
to be used in numerical applications. Since then they have grown into an 
extremely powerful technique for the numerical solution of wide classes of 
nonlinear problems. 

One of the problems to which continuation techniques are applied 
concerns the solution of a nonlinear equation F(x) = in some space 
X. In order to compute a specific solution x* X a possible approach is 
to imbed the equation into a family of equations H (x, t) = 0, < t < I, 
for which there exists a computable solution path x = x(t), 0<t<l, in 
X beginning at a known point z(0) and ending at the desired solution 
z(l) = x*. In other words, one considers a family of smoothly changing 
problems, the final problem being the original problem in question and the 
initial problem being one whose solution is easily determined. Use of the 
continuation, therefore, requires an ability to solve a sequence of problems 
when the solutions to nearby problems are known. 

A related, but conceptually distinct, problem arises hi many applica- 
tions in science and engineering where the equation under consideration 
always includes a number of physically relevant, intrinsic parameters. In 
other words, the equation has the generic form F(x,p) = 0, where x belongs 
to some state space X and p varies in a parameter space P. In this setting 
it is rarely of interest to determine only a few specific solutions. Instead, 
the requirement is to follow paths on the solution manifold {(z,p) e X X 
P} F(x, p) = 0} in the space XxP of all state and parameter variables and 
thereby to detect specific features of the manifold that signify a change of 
behavior of the system under study. 

For the first of these problems, namely the computation of a specific 
solution of some question, two distinct continuation techniques are avail- 
able. The first involves the case when the path of solutions to the family 
is smooth, which in turn allows its representation as a solution of a 
differential equation. The second approach is based on homological rather 
than homotopic concepts and makes use of a numerical form of a result in 
algebraic topology (Speraer's lemma). This approach has been refonnu- 
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lated and is now considered principally in the form of algorithms involving 
piecewise-linear solution paths. Much of the recent research in this area 
concerns these piecewise-linear continuation methods and their application 
to fixed-point problems in economics optimization and game theory. 

Continuation methods for following paths on the solution-manifold of 
parametrized equations developed mainly in structural engineering under 
the name of incremental methods. Applications to other areas, as, for 
example, to computational fluid dynamics and to phase transitions in 
statistical physics have only begun to appear relatively recently. For a 
numerical analysis of a given solution, manifold continuation methods have 
to be considered hi a broader sense as a collection of numerical procedures 
for completing a variety of tasks, including the following: 



1. Follow numerically any curve on the manifold specified by an a 
priori given combination of parameter values with one degree of freedom. 

2. On any such curve determine the exact location of target points 
where a given state variable has a specified value. 

3. K desired, at any such target point switch to the trace of a different 
solution path; 

4. On any such curve identify and compute the initial points where 
stability may be lost; 

5. From any one of the critical points follow a path in the critical 
boundary. In contrast to the case under 1, such paths are typically 
specified by combinations of the parameters with two degrees of freedom 
together with the implicit requirement that all points of the path are in 
the critical boundary. 

6. On any solution path determine the location of bifurcation points 
and the paths intersecting at that point. 



For specific applications additional tasks may arise. For instance, if 
the parametrized equation represents the equilibrium equation of a system 
of differential equations then we may wish to locate Hopf bifurcation points 
on a particular solution path where periodic solutions of the dynamical 
system branch off from the static equilibrium. 

In recent years much research has been devoted to these various tasks, 
but there remain many open questions especially hi connection with the 
more specialized tasks, such as the location of Hopf bifurcation points. 
Moreover, library-quality packages for the basic tasks 1-3 are still under 
development, and the software is by far not so advanced as, say, in the case 
of software for the solution of ordinary differential equations. For specific 
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applications, as, for example, to fluid-dynamic problems or to the case 
of structural problems involving plasticity, creep, and viscoelastic effects, 
the situation is still wide open; and only relatively ad hoc techniques are 
available. 



4.5 OPTIMIZATION METHODS 

Optimization problems occur in all areas of science, engineering, 
economics, and applied statistics. They may involve some least-squares 
approximation of observed data, fitting of parameters occurring in a math- 
ematical model on the basis of experimental observations, optimization of 
the design of an engineering structure, optimal control of an engineering 
system, or modeling of economic or business systems. These are only a 
few examples that do not even begin to cover the numerous types of op- 
timization problems that arise in practice. 

Broadly speaking, in all of these problems a real function, usually 
called the objective function, is to be minimized or maximized over some 
constraint set in a given finite- or infinite-dimensional space. The problems 
differ considerably depending on the mathematical characteristics of the 
objective function and the constraints, the dimension of the space, the 
amount of computable information that is available, and the requirements 
of specific applications. 

A good overview of the research problems concerning computa- 
tional methods for such optimization problems in finite-dimensional spaces 
has been given in the report Program Directions for Computational 
Mathematics, June 1979, prepared for the Applied Mathematical Sciences 
Research Program of the U.S. Department of Energy (edited by R. E. 
Huddleston, Sandia National Laboratories). This report identified the fol- 
lowing research areas: 



1. Analysis and comparison of techniques for dealing with nonlinear 
inequality constraints. 

2. Construction of algorithms for environments where computer 
storage is restricted. 

3. Production of high-quality software and other software-directed 
activities. 

4. Development of algorithms and software for problems in which 
some of the variables are discrete or integer. 
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5. Techniques to assist in finding global (or specified) minima. 

6. Large-scale linear programs with special structure, e.g., staircase 
and block-angular. 

7. Methods for nondifferentiable optimization. 

8. Techniques for selected nonlinear estimation problems, such as 
separable or constrained least squares, and data fitting in special norms. 

9. Investigation of the connection between fixed-point methods and 
alternative methods for solving nonlinear equations. 

10. Further study of selected aspects of linear constrained problems, 
such as special types of constraints, application of conjugate gradient 
techniques, the merits of various active set strategies, and the influence 
or scaling. 



The report does not address infinite-dimensional problems, for ex- 
ample, problems in control theory or the calculus of variations. Here the 
mathematical problems often become formidable, and special techniques 
are needed in most applications. For the computational solution, the prin- 
cipal approach applies finite-dimensional optimization methods to a dis- 
cretized form of the problem. This, in turn, raises the usual questions 
concerning the convergence and quality of the approximations. 



4.6 GRAPH-THEORETICAL METHODS 

Many physical systems, or their mathematical models, involve a number 
of discrete objects that are the producers or users of certain data or 
commodities and that in turn are interconnected by links that can transmit 
these items from one of the objects to another one. The underlying 
mathematical structure is then a graph consisting of a collection of vertices 
connected by arcs or branches. 

The modeling of some classes of physical systems in terms of graphs 
is rather natural, for example, in the case of communication or transpor- 
tation networks where the branches represent, say, power lines, roads, air- 
line routes, or water pipes, and, correspondingly, the vertices are power- 
generation stations, communities, airline terminals, or water reservoirs. In 
other cases the connection between a particular problem and its graph- 
theoretical interpretation is less obvious, for example, in the case of a 
discretization of a particular boundary-value problem on a given grid or 
the representation of a collection of data in a computer. 
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Corresponding to the wide range of applications, the forms of graph 
problems differ widely. For certain network models interest may center 
on the connectivity properties of the graph, that is, on the determination 
of whether a particular commodity can be transported between specific 
vertices. Then questions arise about the maximum possible flow that can 
be accommodated under particular constraints. These problems in turn are 
related to the so-called vulnerability and reliability problems for networks. 
On the other hand, if time enters into consideration then questions of 
waiting time and of best routing may arise. These are only a few of many 
types of such problems. 

From a computational viewpoint these various problems have stimu- 
lated the development of numerous combinatorial and graph-theoretical al- 
gorithms. But there remain many open research questions, especially when 
it comes to the production of general, robust software and the availability 
of algorithms for problems involving large graphs. 

Graph-theoretical formulations are also finding increasing application 
in connection with the numerical solution of problems that do not have an 
inherent graph-structure. One such class of problems concerns the com- 
putational handling of large sparse matrices. Many of the algorithms used 
here perform a sequence of steps that transform the matrix into succes- 
sively simpler matrices of, say, a more nearly upper triangular or diagonal 
type. These transformations achieve their aim by introducing zeros in 
place of originally nonzero matrix elements. But, as an unavoidable by- 
product they also introduce nonzero elements hi places where the original 
matrix elements were zero. Thus, in order to exploit the sparsity struc- 
ture of the matrix one must provide a data structure that either allocates 
from the outset space for all the fill-in during the computation or that 
allows for a dynamical allocation of space for the fill-in when it occurs. In 
both cases, graph-theoretical approaches form the basis for the design of 
desirable algorithms. 

An example of the use of a static data structure is the case when a row 
and column representation can bring the matrix into banded form with 
small bandwidth. This is frequently the case of the linear systems arising 
hi the finite-element discretization of elliptic boundary-value problems, 
and the corresponding bandwidth optimization routines are widely used 
in practice. On the other hand, if a dynamic storage structure is used, then 
special care has to be taken in the design of pivoting strategies for reducing 
fill-in while at the same time maintainmg numerical stability. Various 
sparse matrix packages have been developed for this purpose. They are 
typcially rather large and complex pieces of software. But the field is still 
under active development, and there remain many open research questions. 
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Iii particular, much still needs to be done on the interrelation between 
sparse matrix algorithms and special-purpose computer architectures. 

In recent years graph-theoretical approaches have also found increas- 
ing application in the design of algorithms for the numerical solution of 
problems of mathematical physics, in particular, of fluid dynamics. For ex- 
ample, it has been observed that certain natural finite-difference discretiza- 
tions of the equations of viscous, incompressible flow admit interpretations 
as systems defining flows on certain associated networks. Typically, the 
network nodes correspond to the idealized control volumes represented by 
the mesh points of the finite-difference equations, whereas the network 
arcs correspond to the paths on which one may identify the discrete finite- 
difference mass fluxes. Such observations can lead to remarkable savings 
in computing costs because they open the way to a priori transformations 
of the original (discretized) system of equations to completely equivalent 
systems in substantially fewer variables. For flow problems in two and 
three space dimensions the reduction factors are nominally 1/3 and 1/2, 
respectively. In graph-theoretical terms this approach corresponds to the 
transformation of the Maxwell-node equations to the Maxwell-mesh equa- 
tions, long used in electrical-circuit problems. 

The key requirement hi the application of network techniques to fluid- 
flow problems is that the discretized momentum and continuity equations 
be interpretable as "Ohm's laws" and "Kirchhoff-node laws" on an as- 
sociated network. This requirement permits great generality in the actual 
form of these laws and easily accommodates both implicit time-dependent 
as well as nonlinear steady-state difference forms of the Navier-Stokes 
equations. 

It is even possible to extend the network approach to compressible flow 
problems. Here the new idea appears to be the introduction of pseudo- 
flows into the node law and the identification of corresponding pseudo- 
characteristics to augment the Ohm's laws obtained from the momentum 
equations. In the nonstationary case, this produces numerical schemes 
whose stability is independent of local acoustic velocities. 

The entire area is still under active research development. For ex- 
ample, there are efforts to bring multispecies flow problems under the 
purview of network theory. This includes, as an important special case, 
two-phase flow problems that are a central concern in the design of reliable 
energy-generating systems and the safety analysis of nuclear power plants. 
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4.7 SPLITTING METHODS AND DEFECT CORRECTIONS 

Splitting is a means to separate a large intractable problem into a sequence 
of smaller, or at least more easily solved, problems. These methods are 
invoked to reduce significantly the computational effort (time and memory) 
needed to solve a problem compared with solving it directly. Often a 
form of splitting is required to make solution practical, and some cases 
even tractable. Various forms of splitting are common in engineering and 
scientific applications, even if they are not always recognized as such. 

Defect correction is a widely used, if often unlabeled technique of 
splitting. It presumes that one wants to solve a given hard problem, that 
one has in hand a guess at its solution, and that one also has a nearby 
problem that can be solved easily. It corrects the guess by solving the 
easier problem with special data computed from the guess. 

To illustrate the method, suppose that after discretization it is neces- 
sary to solve a finite-dimensional system of equations of the form 



- 

where A is a nonlinear operator, 6 is a known vector, and v is the solution. 
Often v is difficult to obtain directly, but the residual error 

r A(w) b 

for an approximate solution w is easy to evaluate. If there is a related 
system 



that approximates the original system and that is easier to solve, the defect 
correction algorithm may be appropriate. Given an estimate v n near a 
root v n+ i of the original system, we can expand the original equation in 
a Taylor series to get 

AK) - 6 + P(v n+1 ) - P(v n ] - ( J P - J A )(v n+l - v n ] + 0(e 2 ) 

where e = v n+ i v n . The defect correction iteration is any approximation 
to the above equation. The simplest such iteration is 



This iteration will converge if v n and J P (the Jacobian of P) are near 
enough to v n+i and J At respectively. 
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One of the most common splitting algorithms in engineering applica- 
tions is the alternating direction implicit (ADI) method. Here, the model 
problem is a partial differential equation system having two or more in- 
dependent spatial coordinates. Its direct solution requires too much com- 
puter storage and time to be tractable. Using an ADI splitting, the prob- 
lem is reduced to a manageable sequence of one-dimensional problems. 
In terms of the defect corrections, the system P is the sequence of one- 
dimensional problems that is much easier to solve, yet that closely ap- 
proximates the original system. 

Other forms of splitting also arise. For example, in certain systems of 
partial differential equations some of the equations are weakly coupled to 
the others. In these cases, solving the equations sequentially (rather than 
coupled) can result in significant savings. Similarly, in some combustion 
problems, considerable savings are realized through operator splitting al- 
gorithms in which the chemical rate terms are handled separately from the 
species transport terms. These methods are equivalent to matrix splittings 
of the Jacobian of the system. 

Analysis of splitting methods is important because the methods often 
do not converge. We must be concerned about accuracy and convergence 
of the factored system. Even though each iteration may be fast, we 
must have some idea about the overall convergence rates and about any 
degradation in accuracy resulting from the split. Analysis will also likely 
lead to the identification of matrix properties that suggest a beneficial 
splitting that would not be apparent from physical reasoning. Splitting 
can be considered as either a splitting of the original equations or as an 
approximate factorization of the iteration matrix. The former is the most 
intuitive, and the one in which physical insight is valuable. However, the 
latter is probably the one more amenable to analysis. 



4.8 MONTE CARLO TECHNIQUES 

Mathematical solution methods can be broken into two broad classes, 
Monte Carlo methods and deterministic methods, depending on whether 
chance events are included in the method. All scientists are familiar 
with deterministic methods but most have little familiarity with Monte 
Carlo methods. Deterministic methods do not depend on chance, and 
calculations performed using the same input data and the same computer 
code will provide exactly the same results. Monte Carlo calculations using 
the same input data and the same computer code will provide different 
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results, depending on what chance events occur. In this section we discuss 
what Monte Carlo methods are and how they might be improved. 

The Monte Carlo method estimates averages on a probability modd. 
Any quantity that can be expressed as an average value can be calculated 
by Monte Carlo techniques. Sometimes a probability model is obvious 
from the problem itself. For example, the probability that tossing three 
fair coins and obtaining a desired outcome [for example two heads (H) and 
one tail (T)] is easily estimated by Monte Carlo methods. The probability 
model consists of assigning a probability of ^ to an H and a probability 
of to a T on each toss and assigning a score of 1 to a desired outcome 
(HHT, HTH, THH) and a score of to any other outcome (HHH, TTT, 
THT, HTT, TTH). The computer has a random number generator that 
generates random numbers uniformly on the unit interval (0,1). A uniform 
distribution means that any random number is equally likely to have any 
value between and 1. Thus a com toss can be simulated on a computer 
by: 

H occurs if < 5 
T occurs if > f 

To toss three coins the computer uses three random numbers &, &, 
and 3 . A typical set of tosses might be i = 0.7 (T), & = 0.1 (H), and 
3 = 0.4 (H), scoring 1. The probability of obtaining two heads and one 
tail is (approximately) the average score (|) after many sets of computer 
tosses. 

Sometimes a probability model is not immediately apparent, but after 
a little thought the desired calculation can often be expressed as the 
estimation of an average value. For example, the integral 



can be thought of as the average value of f(x) on the interval (a, b). To 
see this, note that by definition the average value of f(x) with respect to 
a probability density p(x] is 



f b 
</>=/ f(x)p(x)dx. 

Ja 



Thus J is the average value of f(x) with respect to the probability density 
p(x} = l/(b a). Here p(x) is a uniform probability density on (a, b), and 
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x values can be samples uniformly on (a, 6) by setting 



whereupon 

i N 



where I is the estimate of I. . , , . , 

There are many applications for the Monte Carlo technique, but the 
application to neutron and photon transport has probably consumed more 
computer time than all other applications combined. Because ol tms, ana 
because neutron and photon transport problems have a natural pro baouiy 
model, the remaining discussion will pertain to neutron or photon ^trans- 
port problems. Furthermore, because neutron and photon transport are 
similar, the remaining discussion will refer to particle transport. 

Particle flow in nature is a random process. One cannot say exactly 
what will happen to an individual particle. One can only say what the 
probabilities are that various events occur. For example, a particle ttavei- 
ing through matter has a certain probability to collide per unit s , 
the actual distance between collisions is unknown, but the probab iihty ol 
traveling a distance I without collisions is known. Smulariy the ^nucuae 
a particle will collide with is not known, but the probability of colliding 

with the nuclide is known. ^,Wo-mQ ns 

The simplest Monte Carlo model for particle transport P^lems uj 
the natural probabilities that various events occur for the probabib 
model in essentially the same way as the coin toss example. That 
particles are followed from event to event in a ComputerLand tne 
is always selected (using the random number generator) trom a 
possible next events according to the natural event ****** 
type of model is called the analog Monte Carlo model because it is 
analogous to the naturaUy occurring ^transport 



particles contribute, and the statistical uncertainty m the answer 



models for particle transport that will estimate the same average value 
as the analog Monte Carlo model. These other techniques are known 
as nonanalog Monte Carlo models, and they are useful because although 
the average value remains unchanged, the variance (uncertainty) of the 
estimate can often be made much smaller than the variance for the analog 
estimate. Practically this means that problems that would be impossible 
to solve in days of computer time can be solved in minutes of computer 
time. 

A nonanalog Monte Carlo model attempts to follow "interesting" par- 
ticles more often than noninteresting particles. An "interesting" particle, 
by definition, is a particle that contributes a large amount to the quantity 
(or quantities) that needs to be estimated. There are many nonanalog 
techniques, and they all are meant to increase the odds that a particle 
scores (contributes). To ensure that the average score is the same in the 
nonanalog game as in the analog game, the score is modified to remove 
the effect of biasing (changing) the natural odds. Thus, if a particle is 
artificially made q times as likely to execute a given random walk, then 
the particle's score is weighted by (that is, multiplied by) 1/q. The average 
score is thus preserved because the average score is the sum, over all ran- 
dom walks, of the probability of a random walk and the score due to that 
random walk. The trick in obtaining low- variance solutions is to choose g's 
such that all random walks contribute the same score, in fact the average 
score. This then would be a zero-variance solution. 

It is always possible for any (linear) particle transport problem to 
select q's for each random walk such that every particle's score is the 
average score; that is, a zero-variance solution. The hooker is, of course, 
that the perfect q'& are not known, thus a zero-variance solution is not 
practical. However, people have often been successful enough in guessing 
good g's, that is biasing the odds, to improve the calculational efficiency 
by several orders of magnitude. This is typically done iteratively with a 
person making several short Monte Carlo computer runs (calculations) and 
using information gained on each run to better guess how to bias the next 
run. When the person's guesses no longer improve the calculation, a long 
run is done with the optimal biasing learned in the short runs. 

Can the computer learn to optimize the biasing? The computer is, 
after all, supplying the information that the person uses to learn. In 
the past decade a number of computer learning techniques have been 
tested. Thus far it has proven impossible for a computer to take an 
arbitrary transport problem and, without human aid, optimize the biasing. 
However, two things are worthy of note. First, computer learning aided 
by human judgment appears to be substantially better in many cases than 
human learning alone. This typically works by having the computer inform 
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the person how the computer would bias the subsequent run and having 
the person selectively alter the computer's suggestions according to the 
person's intuition. Second, the amount of human judgment required is 
decreasing. The day may come when the computer needs no human aid. 
Once human aid is no longer needed, the computer can learn to 
adjust its own biasing particle, as the calculation proceeds. An interesting 
implication of an adaptive Monte Carlo technique is that the common 
central limit theorem of statistics that would constrain the accuracy of 
the calculation to decrease as the square root of the number of particles 
followed no longer applies. The common central limit, theorem applies 
only when each particle's random walk is independent of all others and 
the sampling process is identical for all particles. Consequently, the rate 
of convergence may be more rapid. Although some limited statistical 
theories exist for dependent random variables, it has not been investigated 
whether these theories can profitably be applied to Monte Carlo particle 
transport problems. Thus, it is uncertain how fast an adaptive Monte. 
Carlo calculation is converging or indeed what the maximum convergence 
rate for a good learning process might be. Empirical evidence shows that 
the convergence can be substantially faster than the square root of the 
number of particles. In light of this, Monte Carlo methods can be expected 
to become more competitive with deterministic calculations. 



4.9 PROBLEM-DEPENDENT METHODS 

A variety of adaptive methods have a common goal: to overcome the 
computer size and grid resolution limitations, which are especially severe 
in singular problems, by use of computational elements that fit or model 
the singularity more directly. In this way, it may be possible to incorporate 
into the solution algorithm considerable analytic information beyond that 
provided by the equations themselves. This theme occurs in many aspects 
of numerical analysis. In the finite-element method, one may choose 
elements that include any singularities in the solution being computed (see 
Section 3.3). Grid adaptivity (see Section 4.2) is also problem dependent. 
Here we discuss examples of somewhat more special methods. Of necessity 
they apply only to classes of problems that possess related singularities. 
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4.9.1 The Riemann Problem and Nonlinear Wave Modes 

The Riemann problem for a nonlinear hyperbolic system of conservation 
laws is the Cauchy problem in one space dimension for data that are 
everywhere constant except for a single jump discontinuity at the origin. 
The solution of the Riemann problem provides a resolution of this discon- 
tinuity into the nonlinear modes, or waves, which propagate coherently 
in time. This idealized problem can be thought of as an approximate 
description of a higher-dimensional flow field especially in the neighbor- 
hood of a discontinuity hypersurface. This point of view has led to a num- 
ber of improved finite-difference schemes, which attempt to determine the 
various nonlinear wave modes at each point in space and time and to ad- 
just the differencing of the differential equations to take advantage of this 
knowledge. This adaptivity is to the structure of state space, in contrast 
to the coordinate space adaptivity discussed in Section 4.2. 



4.9.2 Front Tracking 

Front tracking is a combination of adaptive grid methods with the use of 
Riemann problems. The method is adapted to problems that have impor- 
tant singularity hypersurfaces (lines in two space dimensions, surfaces in 
three space dimensions), such as shock waves and contact discontinuities. 
In one version of this method, there is an overall time-dependent coordinate 
transformation to map the discontinuity into a fixed grid line in some set 
of computational coordinates. In another version, the discontinuity is rep- 
resented by a moving lower-dimensional grid that follows ("tracks") the 
dynamical motion of the discontinuity. The motion of the discontinuity 
and of the moving grid points that track it are governed by solutions of 
Riemann problems, or equivalently by the method of characteristics. 



4.9.3 Vortex Method 

This method introduces small lines or surfaces of vorticity into a fluid flow. 
The method is adapted to the study of shear-layer discontinuities, bound- 
ary layers, boundary-layer separation, and turbulence. The equations of 
motion of an ideal fluid yield simple equations for the motion of a collec- 
tion of vortices imbedded in the fluid. In fact, the vortices move passively 
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with the fluid, and their mutual interaction is described by a Hamiltonian 
system of ordinary differential equations, with Coulomb type interaction 
energy. La the case of the Navier-Stokes equations, the vortex motion also 
contains a diffusion term. 

These methods have been applied successfully to the problem of tur- 
bulent flame propagation (see Section 2.2). In this problem, the turbulent 
mixing is a primary factor in determining the flame speed. The turbulence 
comes from the boundary layer and in the boundary layer is calculated by 
a vortex method. Related methods have been developed under a variety of 
names such as boundary integral methods and Green's function methods. 



4.9.4 Scale Invariance and the Renormalization Group 

Scale transformations are the transformations x * sx, acting on space or 
on space and tune. A function u is homogeneous of degree a if 

u(sx, ay, sz) = s a u(x, y, z) 

and scale invariant if a = 0. Many problems have solutions that are scale 
invariant or approximately scale invariant over some parameter range. 
Such solutions are called similarity solutions. Using the scale invariance, 
one of the independent variables can be eliminated, making the solution 
more elementary to compute. 

However, scale invariance can also indicate the occurrence of com- 
plex phenomena. Specifically, any singularity that occurs in a scale- 
invariant problem must be repeated in all length scales (for which the 
scale invariance holds). Mathematically, Cantor sets, snowflake curves, 
and fractals are examples of such phenomena. In statistical physics, criti- 
cal phenomena in the equation of state is a scale-invariant problem. One 
general picture of turbulence holds that scale invariance (vortices on a 
large range of length scales) describe the inertia! range, or energy transport 
range, of turbulence. 

To implement these ideas in a computational algorithm, one integrates 
over a given set of length scales in a standard manner. The result of this 
computation is taken as data for a new calculation over a new set of length 
scales (with the original degrees of freedom eliminated from the prob- 
lem). This process is iterated and if convergent gives the scale exponent 
a, Sample numerical calculations of this type were discussed in Section 
2.1 in connection with turbulence. The method is well established for per- 
turbative calculations (with a small parameter) in critical phenomena in 
statistical physics. 
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FIGURE 4.1 The evolution and merger of isolated vortex structures as 
predicted by contour dynamical techniques. [From E. A. Overman n and 
N. J. Zabusky, Phys, Fluids 25, 1297-1305 (1982).] 
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Since the renormalization group methods are novel, it is worth men- 
tioning that the mathematical foundations of this method have been estab- 
lished recently in several cases including examples of hierarchical models 
in statistical mechanics, interval maps, and renormalized quantum fields. 



4.9.5 Contour Dynamical Methods 

Contour dynamical methods are being applied to a variety of incompress- 
ible flows in two dimensions. These generalizations of the "waterbag" 
method provide simplified models for following the evolution of contours 
separating regions of piecewise-constant (pc) density that are the sources 
of the flow. The flow is the result of the self and mutual interaction of con- 
tours that evolve, mainly by area-preserving maps. These methods have 
been applied to the Euler equations, where pc finite-area-vortex regions 
and/or vortex sheets at density interfaces are sources of the flow; and 
the equations for a weakly ionized and strongly magnetized ionospheric 
plasma cloud, where pc ion-density regions are sources of the flow. For the 
former, a large class of steady-state translating and rotating solutions with 
pc vorticity ("V-states") have been found. Figure 4.1 shows the merger 
and breaking for a perturbed corotating V-state (two pc finite-area vortex 
regions with the same circulation), a familiar process in free-shear layer 
experiments. Notice how the two regions merge to form one region (by 
snipping out the common boundary at t = 10) and then stabilize by eject- 
ing vorticity in thin filaments. With these methods it is possible to obtain 
detailed information about the regions because the dimension is reduced 
by one. The curvature of the contour provides a predictive signature of 
the evolving small-scale structure, e.g., the roll up of vorticity filaments, 
etc. 



4.10 COMPUTER SOFTWARE AND HARDWARE 

As previously discussed, large-scale scientific calculations that tax the 
resources of the most powerful computers will continue to be essential to 
modern research and development efforts. To obtain long-term reliability 
and stability of future applications codes, implementing and testing of 
high-level numerical software should be coordinated. This will require a 
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strong research and development effort with cooperation supported among 
applications programmers, the theoretical numerical methods researchers, 
and the computer-science software developers. 

Repetition of expensive, error-prone, and time-consuming coding of 
commonly used methods should be avoided. Much of the current scientific 
software now being developed is redundant. If the common elements of the 
existing codes were available as modules, future applications programmers 
could use these routines and eliminate much of their efforts. New software 
is most readily accepted if it is compatible with existing techniques and 
simple enough so that potential users can observe tangibly better results 
in a trial run than those existing methods can produce. If such high- 
level routines were available, they could perform many of the common 
procedures found in applications codes, including grid generation, rezon- 
ing, numerical interpolation, differentiation, and integration; they could 
approximate differential boundary conditions and solve large, sparse non- 
linear systems of equations. 

An important goal is the machine independence of applications pro- 
grams. This is a difficult task because methods tailored specifically for a 
particular machine architecture will probably become more the norm than 
the exception. We can, by keeping machine-dependent codes in libraries of 
high-level software with standard user interfaces, strive to keep the user's 
scientific applications codes portable. The underlying library routines will 
be, necessarily, less portable because the architecture of the new machines 
will certainly be different from that of today's supercomputers. To uti- 
lize the inherent powers of the new machines we will have to re-examine 
traditional methods and identify the better ones for a particular machine 
architecture. 

The continuing revolution in microelectronics is having a profound 
impact on scientific computing. Indeed, it is likely to change our concept 
dramatically of what scientific computations can and cannot be effectively 
performed. Most certainly the impact of this revolution will be much 
greater than, say, the impact that floating-point hardware had. Moreover, 
while the costs of individual tasks will be greatly reduced, the domain of 
scientific computations will be greatly expanded and "frontier computa- 
tions" will continue to be expensive. 

These changes are being brought about by a number of factors: 
Individual components are becoming increasingly faster and smaller. Very- 
large-scale integration of such components is not only reducing the size of 
the packaged systems but also providing opportunities for customized, in- 
formation processors. Also, the decrease in the cost and size of computer 
memory implies that we can look for much larger memory systems. This 
will obviate many of the existing difficulties with secondary storage. 
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Despite the fact that components are becoming faster, the limits of 
raw machine speed are not visible, and further gains will eventually have 
to be made by using clever architectures and algorithms. Some kind of 
parallelism seems to be unavoidable. The programming issues involved in 
using parallel machines are still not satisfactorily resolved. The automatic 
detection of parallelism and the resulting scheduling of multiple resources 
are important open problems. 

Architectures, such as systolic arrays, based on specific subtasks can 
increase the performance of systems involving these subtasks by several 
orders of magnitude and clearly have a bright future. Algorithms contain- 
ing compute-intensive subtasks that can be vectorized in this fashion have 
a bright future. Because these architectures are hi general regular, the 
algorithms that can be vectorized for such machines tend to be regular, 
i.e., simple, nonadaptive, uniform-grid, low-order algorithms. It is clear 
that there are nicely behaved problems for which these regular algorithms 
on specialized machines will require significantly less time than algorithms 
requiring fewer operations on serial machines. It is also clear that no 
matter how fast the specialized machines are, there are problems that 
are sufficiently difficult that more sophisticated algorithms are needed for 
more general-purpose computers. 

In order to bring about these advances in architecture, it is necessary 
to involve practitioners of scientific computation in the design process. 
Luckily, modem design automation tools should make it possible for in- 
terdisciplinary design teams to successfully synthesize innovative special- 
purpose systems, and automated fabrication facilities should make it pos- 
sible for such systems to be built, debugged, and used. 

These advances in machine architecture should also have a dramatic 
effect on the design and analysis of algorithms for scientific computing. 
Traditionally, such analysis has been based on (asymptotic) estimates of 
the number of arithmetic operations. However, with these new architec- 
tures it is quite likely that the running time of an algorithm will be more 
dependent on the movement of data than on the number of arithmetic 
operations. Thus, we need to develop new analytic models of complexity of 
scientific algorithms so that such models give us useful information about 
the relative performance of algorithms. 
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